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Abstract 

The  Unified  Flow  Solver,  a  hybrid  continuum-rarefied  code,  is  used  to  investi¬ 
gate  the  internal  structure  of  a  normal  shock  wave  for  a  Mach  range  of  1.55  to  9.0 
for  Argon,  and  1.53  to  3.8  for  diatomic  Nitrogen.  Reciprocal  shock  thickness,  den¬ 
sity,  temperature,  heat  flux,  and  the  velocity  distribution  function  are  calculated  for 
a  one- dimensional  shock  wave  and  compared  with  experimental  data  from  Alsmeyer 
and  DSMC  results  from  Bird.  Using  the  Euler,  Navier-Stokes,  BGK  model,  and 
Three- Temperature  BGK  model  schemes,  results  from  UFS  compare  well  with  exper¬ 
iment  and  DSMC.  The  Euler  scheme  shows  atypical  results,  possibly  resulting  from 
modifications  made  to  include  internal  energies. 

An  entropy  spot  is  introduced  into  a  two-dimensional  domain  to  investigate 
entropy-shock  interactions  over  a  range  of  Knudsen  numbers  (Kn  =  0.01,  0.1,  and 
1.0)  for  Mach  2.0  in  Argon.  Previous  work  on  entropy-shock  interactions  has  only 
been  performed  using  an  Euler  scheme.  Here,  results  are  presented  in  Argon  using 
coupled  BGK  and  Navier-Stokes  solvers.  Density,  pressure,  and  temperature  profiles, 
as  well  as  the  profiles  of  their  gradients,  are  reported  at  certain  times  after  the  entropy 
spot  convects  through  the  shock. 
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An  Investigation  of  Shock  Wave  Physics  via  Hybrid 
CFD-BGK  Solution  Methods  for  Nonequilibrium  Flows 

I.  Introduction 

Hypersonics  is  a  current  field  of  research  and  interest  in  the  Department  of 
Defense  (DoD),  the  United  States  Air  Force  (USAF),  and  around  the  world.1  For 
example,  in  2006  the  Hypersonic  International  Flight  Research  Experimentation  pro¬ 
gram  (HiFire)  began  with  $54  million  of  funding  and  an  agreement  between  USAF 
and  the  Australian  Department  of  Defense.  HiFire  will  continue  through  2012,  pro¬ 
viding  basic  and  applied  research  to  better  understand  hypersonic  flows,  as  well  as 
to  generate  experimental  data.  The  program  will  include  up  to  10  flight  experiments 
(with  payload)  at  realistic  hypersonic  flight  conditions.  Seen  as  a  game-changing  ca¬ 
pability,  Douglas  Dolvin  of  the  Air  Vehicles  Directorate  of  the  Air  Force  Research 
Laboratory  (AFRL)  stated, 

We  envision  air-breathing  powered  hypersonic  cruise  missiles  in  the  near- 
term,  which  are  able  to  deliver  prompt,  precision  strike  of  time  critical 
targets  from  safe,  standoff  distances.  In  the  far-term,  these  air-breathing 
hypersonic  vehicles  may  enable  operationally  responsive  space  access.  [1] 

In  order  to  understand  the  difficulties  associated  with  hypersonics  research,  it 
is  important  to  understand  more  about  hypersonic  flows.  A  flow  usually  defined 
as  hypersonic  if  the  Mach  number  is  greater  than  or  equal  to  five  (  M  >  5  ), 

although  Anderson  has  extended  the  definition  to  include  any  flow  where  certain 
physical  phenomena  become  important  [4], 2  Civilian  transport  aircraft,  such  as  a 
Boeing  747,  operate  in  subsonic  speeds  (  M  <  1  )  and  up  to  a  certain  altitude.  Due 
to  design  concerns,  such  aircraft  do  not  break  the  sound  barrier,  and  in  practice  stay 
below  transonic  regimes  (  0.8  <  M  <  1.2  ),  where  some  flow  over  the  aircraft  is 

1Lists  of  symbols  and  abbreviations  are  included  in  the  appendices. 

2Such  as  a  thin  shock  layer,  an  entropy  layer,  strong  viscous  interaction,  high  temperatures,  and 
low  densities. 
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supersonic  (  M  >  1  )  and  some  of  the  flow  is  still  subsonic.  On  the  other  hand, 
fighter  jets,  such  as  the  F-22  Raptor,  are  designed  to  fly  at  supersonic  speeds,  and 
regularly  do  so  when  flying  over  oceans.  At  supersonic  speeds,  shock  waves  form  due 
to  the  compressibility  effects  of  the  air.  At  such  high  velocities,  the  kinetic  energy 
of  the  flow  can  be  dissipated  in  the  form  of  heat  due  to  viscous  effects,  thereby 
significantly  raising  the  temperature  of  the  aircraft  surface.  With  increasing  velocity 
(and  Mach  number),  the  kinetic  energy  of  the  air  increases,  making  dissipative  heat 
transfer  more  of  a  limiting  factor,  to  the  point  where  protruding  surfaces  (such  as 
control  surfaces  required  for  steering)  are  extremely  susceptible  to  thermal  failure. 
As  a  result,  hypersonic  vehicles  tend  to  have  a  more  streamlined  appearance,  such  as 
the  space  shuttle,  which  is  a  atmospheric  re-entry  vehicle.  Besides  thermal  concerns, 
hypersonic  vehicles  require  sophisticated  propulsion  systems  and  must  operate  at  very 
high  altitudes  in  order  to  avoid  the  higher  air  densities  that  exist  near  sea  level.  A 
higher  air  density  corresponds  to  more  molecules  per  unit  volume,  leading  to  increased 
drag  during  flight.3  Therefore,  it  becomes  prohibitive  to  achieve  hypersonic  flight 
except  in  rarefied  regions.  It  now  becomes  more  obvious  why  hypersonics  research  is 
very  difficult  and  expensive  to  perform. 

The  most  common  way  to  avoid  the  high  cost  of  experimental  hypersonics  work 
is  to  model  such  flows  using  computational  codes.  Computational  Fluid  Dynamics 
(CFD)  codes  allow  continuum  flows  to  be  analyzed  using  the  Euler  and  Navier-Stokes 
(NS)  equations.  As  will  be  explained  in  later  sections,  the  Euler  and  NS  equations 
are  not  capable  of  accurately  describing  nonequilibrium  flows.  In  fact,  the  Euler 
equations  assume  the  flow  is  at  equilibrium,  while  the  NS  equations  assume  only 
small  deviations  from  equilibrium.  Also,  the  Euler  and  NS  equations  are  derived 
assuming  that  the  flow  is  a  continuous  mass  of  fluid  (continuum),  which  can  no  longer 
be  assumed  at  low  densities  (rarefied).  Therefore,  CFD  codes  breakdown  in  both 
nonequilibrium  and  rarefied  flows.  Kinetic  methods,  however,  make  no  assumptions 

3The  term  molecule  here  refers  to  both  monatomic  and  polyatomic  particles.  This  usage  will  be 
followed  throughout  the  rest  of  this  work. 
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about  density,  but  rather  about  how  the  separate  molecules  interact.  These  methods 
are  commonly  derived  from  kinetic  theory,  starting  at  the  molecular  level  with  the 
goal  of  describing  macroscopic  properties  (such  as  density,  pressure,  temperature, 
and  viscosity)  with  molecular  interactions.  Whether  statistically  or  deterministically 
based,  kinetic  methods  usually  describe  the  flow  with  distribution  functions,  such  as 
the  velocity  distribution  function  (VDF).  The  VDF  describes  the  likelihood  that  any 
molecule  in  that  region  will  have  a  given  three-dimensional  velocity.  Alternately,  the 
VDF  can  be  seen  as  the  actual  distribution  of  velocities  over  all  the  molecules  in  a 
given  region.  The  former  statistical  interpretation  can  be  applied  at  any  point ,  while 
the  latter  can  only  be  applied  over  a  region  small  enough  when  compared  with  the 
entire  flow,  but  large  enough  that  it  contains  enough  molecules  that  the  VDF  is  a 
continuous  function.  The  Boltzmann  equation,  as  will  be  seen  later,  is  a  governing 
equation  for  the  VDF,  and  describes  its  evolution  in  time,  physical  space,  and  velocity 
space,  due  to  intermolecular  collisions.  Even  though  the  Boltzmann  equation  has  its 
limitations  (e.g.  only  considers  bimolecular  collisions),  it  has  been  found  to  be  very 
useful,  especially  with  modern  computational  capabilities. 

Many  flows  considered  in  aerodynamics  are  in  a  state  of  equilibrium,  mean¬ 
ing  that  macroscopic  flow  properties  do  not  change  over  time.  However,  there  are 
instances  when  significant  nonequilibrium  effects  come  into  play.  For  example,  con¬ 
sider  a  chemical  process  that  involves  chemical  nonequilibrium.  The  original  chemical 
composition  begins  at  some  initial  state  (usually  at  equilibrium),  then  a  chemical  re¬ 
action  occurs  during  which  the  molecules  are  in  a  state  of  change  (nonequilibrium), 
and  finally  the  chemical  composition  arrives  at  (or  relaxes  to)  a  final  state  (again  at 
equilibrium).  Thus,  nonequilibrium  is  the  process  by  which  a  system  changes  from 
one  state  of  equilibrium  to  another  state  of  equilibrium.  But  there  are  more  types  of 
nonequilibrium  than  just  chemical.  A  flow  may  be  in  chemical  (to  include  dissocia¬ 
tive,  recombinative,  and  ionization),  thermal,  translational,  rotational,  vibrational, 
electronic,  and/or  radiative  nonequilibrium. 
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The  translational  energy  of  a  molecule  is  directly  related  to  its  velocity.  This  en¬ 
ergy  is  transferred  from  one  molecule  to  another  through  a  molecular  collision.  It  only 
takes  a  handful  of  collisions  for  a  molecule’s  translational  energy  to  equilibrate  with 
that  of  its  neighbors,  so  translational  nonequilibrium  is  usually  a  short-lived  process 
(except  for  the  case  of  very  low-density  flows) .  A  measure  of  how  far  apart  molecules 
are  in  space  is  the  mean  free  path  A  (m),  defined  as  the  average  distance  a  molecule 
can  travel  before  it  will  collide  with  another  molecule.  Translational  nonequilibrium 
can  therefore  be  studied  in  flows  with  length  scales  on  the  order  of  the  mean  free 
path,  such  as  shock  waves  or  low-density  (or  rarefied)  flows.4 

A  polyatomic  molecule  has  some  amount  of  rotation  associated  with  it,  leading 
to  its  rotational  energy.  Various  models  are  available  to  describe  how  a  molecule 
rotates,  but  for  this  work  it  is  assumed  that  diatomic  Nitrogen  can  be  modeled  as 
two  Nitrogen  atoms  attached  by  a  rigid  connector.5  Therefore,  N2  has  two  degrees  of 
rotational  freedom,  meaning  that  its  rotational  energy  is  due  to  rotation  about  only 
two  orthogonal  axes  (no  energy  is  associated  with  rotation  about  the  axis  which  passes 
through  both  atoms).  Similar  to  translational  energy,  it  only  takes  on  the  order  of 
10  molecular  collisions  to  reach  rotational  equilibrium.  Most  aerospace  applications 
in  polyatomic  flows  will  involve  rotational  energy  since  rotational  levels  begin  to  be 
excited  at  a  temperature  of  about  2  K.6  It  has  been  observed  that  molecular  rotation 
causes  a  centrifugal  force  that  can  decrease  the  energy  required  for  dissociation.7 

Vibrational  and  electronic  energy,  dissociation,  and  ionization  will  not  be  impor¬ 
tant  in  this  study,  but  are  interesting  nonetheless.  Vibrational  energy  is  a  diatomic  or 
polyatomic  phenomenon.  For  N2,  the  atoms  are  separated  by  some  link,  and  they  tend 
to  vibrate  with  some  energy.  Vibrational  equilibrium  requires  many  more  molecular 
collisions  than  translational  and  rotational  equilibrium  (about  25  trillion  collisions  at 
room  temperature).  However,  vibrational  excitation  for  N2  only  becomes  important 

4Strong  shocks  have  a  thickness  of  only  a  few  mean  free  paths. 

5This  model  is  called  the  rigid  rotator  model. 

6For  JV2,  the  characteristic  temperature  for  rotation  is  2.9  K 

7Dissociation  occurs  when  an  atom  breaks  free  from  the  bonds  that  hold  it  to  the  other  atom(s). 
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at  about  3,390  K.  Electronic  energy  only  contributes  a  small  amount  to  the  overall 
energy  in  practical  applications,  and  is  usually  neglected,  as  is  the  case  here.  Disso¬ 
ciation  and  ionization  are  only  encountered  in  very  high-temperature  flows,  since  in 
N2  they  become  important  at  about  113,000  K  and  181,000  K,  respectively. 

Continuum  flows,  as  described  by  the  NS  equations,  only  experience  small  de¬ 
viations  from  equilibrium.8  Efforts  have  been  made  to  therefore  distinguish  between 
continuum  and  non-continuum  (or  rarefied)  flows,  and  are  in  fact  still  underway  with 
recent  interest  in  hybrid  codes.  One  common  measure  of  whether  a  flow  is  continuum 
or  rarefied  is  the  Knudsen  number  Kn,  defined  as  the  ratio  of  the  mean  free  path  A 
to  a  characteristic  length  of  the  flow  L  (m): 

Kn  =  j  (1) 

A  low  Knudsen  number  corresponds  to  continuum  flow,  while  a  large  value  corresponds 
to  rarefied  flow,  and  a  very  large  value  indicates  free  molecular  flow.  The  choice  of  L 
is  somewhat  arbitrary,  and  really  depends  on  the  problem  at  hand.  For  example,  flow 
over  an  aircraft  wing  might  use  the  mean  aerodynamic  chord  for  L  (small  Kn) ,  while 
the  boundary-layer  flow  over  the  same  wing  could  use  for  L  the  boundary  layer  height 
5gg%  (intermediate  Kn),  and  the  flow  through  the  shock  standing  off  from  the  nose 
would  use  the  shock  thickness  for  L  (large  Kn) .  On  the  other  hand,  A  is  dependent 
mainly  on  the  density  of  the  flow.  The  mean  free  path  would  be  much  larger  at  higher 
altitudes  (high  Kn,  low  density)  than  at  sea  level  (low  Kn,  high  density). 

A  notional  categorization  of  flows  based  on  the  local  Knudsen  number  is  shown 
in  Figure  1.  Continuum  regions  roughly  have  Kn  <  0.01  ,  near-continuum  regions 
span  0.01  <  Kn  <  0.1  (the  continuum  assumptions  begin  to  break  down),  rarefied 
regions  typically  have  values  of  0.1  <  Kn  <  10  (the  continuum  assumptions  are  no 

8Attempts  are  still  being  made  to  extend  the  range  of  applicability  of  continuum  solvers.  For 
example,  Claycomb  studied  how  including  bulk  viscosity  into  a  NS  solver  can  extend  the  effectiveness 
of  a  continuum  solver  into  regimes  with  high  degrees  of  nonequilibrium  [8]. 
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Figure  1:  Notional  local  Knudsen  number  flow  categories. 

longer  valid),  and  flows  with  Kn  >  10  experience  so  few  molecular  collisions  that 
they  are  seen  as  near- collisionless  or  collisionless. 

The  Knudsen  number  is  not  the  only  measure  of  continuum  breakdown.  As  was 
previously  discussed,  continuum  assumptions  are  closely  tied  to  a  state  of  equilibrium, 
or  even  near-equilibrium.  Likewise,  rarefied  regions  display  a  greater  departure  from 
equilibrium  due  to  the  large  time  factors  between  molecular  collisions,  and  therefore 
are  typically  non-equilibrium  flows.  The  study  of  classical  thermodynamics  relates 
the  process  by  which  one  equilibrium  state  changes  to  another  equilibrium  state  with 
a  scalar  parameter,  entropy.  Schrock  and  Carr  investigated  this  further  by  look¬ 
ing  at  how  a  continuum  breakdown  parameter  based  on  entropy  could  be  used  to 
better  characterize  continuum  flow  [7,  34] .  One  drawback  of  their  method  is  that 
the  more  computationally  expensive  rarefied  methods,  such  as  the  Direct  Simulation 
Monte  Carlo  (DSMC)  method,  must  first  be  used  to  calculate  the  entropy  in  order  to 
determine  whether  the  continuum  assumptions  hold  for  that  cell.  Kolobov  et  al.  sug¬ 
gests  two  continuum  breakdown  criteria,  or  switching  parameters  for  the  UFS  hybrid 
code,  based  on  either  density  gradients  or  pressure  and  velocity  gradients  [21].  Using 
the  correct  switching  parameter  is  important  in  UFS,  because  otherwise  non-positive 
VDFs  may  be  obtained,  which  is  physically  impossible  [21]. 
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Besides  determining  continuum  and  rarefied  (or  kinetic)  domains,  hybrid  codes 
must  also  choose  an  appropriate  coupling  method.  Three  classes  of  coupling  have 
been  considered  to  date  [21].  The  first  decomposes  the  physical  domain  into  kinetic 
and  continuum  sub-domains.  The  second  instead  decomposes  the  velocity  domain 
and  considers  separately  fast  and  slow  particles.  The  third  calculates  the  VDF  in 
all  the  cells  and  then  uses  that  information  to  compute  transport  properties  for  the 
continuum  equations.  UFS  takes  the  last  approach. 

There  is  a  trade-off  between  efficiency  and  accuracy  in  computational  modeling. 
Consider  as  an  example  of  extreme  accuracy  the  flow  over  an  airfoil  at  sea  level.  Here 
the  ambient  density  is  relatively  high,  which  means  that  molecules  are  spaced  closely 
together.9  Using  classical  mechanics,  one  could  conceivably  follow  the  velocity  and 
position  (in  spite  of  the  Heisenberg  Uncertainty  Principle)  of  every  molecule.  Every 
time  a  molecule  collides  with  another  (assuming  a  billiard-ball,  or  hard-sphere  model), 
the  deflection  angles  and  momentum  changes  would  be  calculated.  Of  course,  air  is 
composed  of  many  different  species  (e.g.  02,  N2,  etc.),  which  would  only  complicate 
things.  The  flow  of  such  a  large  amount  of  molecules  over  even  a  simple  airfoil  would 
be  prohibitively  expensive  to  compute.  Besides,  how  would  one  know  the  exact  start¬ 
ing  positions  and  velocities  of  every  particle?  In  addition,  the  hard-sphere  model  of  a 
molecule  ignores  all  internal  energies,  only  accounting  for  an  effective  collision  diam¬ 
eter,  molecular  mass,  and  translational  energy.  Instead  of  the  previous  deterministic 
approach,  one  might  take  a  statistical  sampling  of  an  aggregate  of  molecules,  and 
approximate  it  as  one  simulated  molecule  (as  is  the  case  with  DSMC).  However,  by 
making  the  problem  simpler  and  more  manageable  with  simulated  molecules,  one  has 
also  lost  a  measure  of  accuracy.  The  main  point  that  must  be  determined  for  any 
given  computational  model  is  how  well  it  balances  this  trade-off. 

In  an  effort  to  more  efficiently  and/or  accurately  describe  rarefied  flows,  the 
Air  Vehicles  Directorate  of  AFRL  recently  contracted  the  CFD  Research  Corporation 

9It  is  estimated  that  at  sea  level,  there  are  2.5  x  1025  molecules  per  m3  of  air,  and  A  is  on  the 
order  of  0.1  pm. 
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(CFDRC)  to  create  a  kinetic-based  hybrid  code  called  the  Unified  Flow  Solver  (UFS) 
[20-22],  UFS  is  meant  to  be  more  efficient  than  a  DSMC  code  for  near-continuum 
and  transition  flows,  and  more  accurate  than  Euler  and  NS  solvers  for  transition  and 
rarefied  flows.  A  hybrid  code  is  one  that  can  solve  rarefied  and  continuum  flows,  as 
well  as  transition  flows  (flows  that  have  both  rarefied  and  continuum  regions  such  as 
the  hypersonic  flow  over  a  cylinder).  UFS  uses  the  Boltzmann  equation  to  solve  for 
the  VDF  in  each  cell,  and  then  calculates  the  macroscopic  flow  quantities  from  the 
VDF.  Hybrid  solvers  continue  to  be  an  active  research  area,  with  more  developments 
expected  in  the  future. 

This  work  has  two  parts.  The  first  studies  one-dimensional  shock  waves  in 
Argon  and  Nitrogen  (A2)  using  UFS  and  comparing  the  results  to  the  experimental 
data  from  Alsmeyer  and  numerical  simulations  using  Bird’s  DSMC  method  [3].  Shock 
structure  provides  a  useful  test  case  for  any  new  code,  such  as  UFS,  since  it  involves 
translational,  rotational,  and  vibrational  molecular  nonequilibrium. 

The  second  part  of  this  study  examines  how  an  entropy  spot  (i.e.  a  tempera¬ 
ture  or  a  density  spot)  affects  a  two-dimensional  shock  in  Argon  using  one  of  UFS’s 
kinetic  solvers  for  Kn  =  0.01  ,0.1,  and  1.0  at  Mach  2.0.  There  are  three  basic 
modes  of  disturbance  in  a  gas  (acoustic  or  pressure,  vortical,  and  entropy),  and  any 
disturbance  can  be  decomposed  into  a  linear  combination  of  these  three  modes  [12], 
Initially,  studies  were  focused  on  noise  generation  (acoustic  effects)  due  to  turbulence 
traveling  through  a  shock,  using  a  small  vortex  perturbation  in  the  flow  to  simulate 
localized  turbulence  [5,15,25,26,28,31,41],  More  recently,  articles  by  Duck  et  al. 
have  extended  the  freestream  disturbance  analysis  to  include  acoustic,  vorticity,  and 
entropy  waves,  showing  analytically  that  when  any  one  of  these  disturbances  travels 
through  an  oblique  shock,  all  three  modes  are  created  downstream  [9,10].  As  dis¬ 
cussed  by  Fabre  et  al.,  a  cylindrical  entropy  spot  can  be  decomposed  using  Fourier 
synthesis  into  plane  entropy  waves  arranged  with  different  orientations  [11],  Thus, 
an  entropy  spot  traveling  through  a  normal  shock  will  excite  acoustic,  vortical,  and 
entropy  disturbances  downstream  of  the  shock. 
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Hussaini  and  Erlebacher  provide  multiple  reasons  why  entropy-shock  interac¬ 
tions  are  important  in  aerospace  applications,  such  as  turbulence  amplification  in 
shock-turbulent  boundary  layer  interactions,  density  fluctuations  in  supersonic  wake- 
shock  interactions,  noise  generation  in  hot  rocket  exhausts  with  oblique  shock  waves, 
and  enhanced  mixing  caused  by  shock  interactions  in  the  combustor  of  a  scramjet  en¬ 
gine  with  hot  and  cold  flows  (oxidant  and  fuel)  [16].  Numerical  studies  have  been  per¬ 
formed  on  entropy-shock  interactions  using  two-dimensional  Euler  schemes  [11-13,16]. 
However,  no  results  have  yet  been  presented  with  either  the  NS  equations  or  with  a 
kinetic  solver.  This  work  hopes  to  begin  such  a  study. 
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II.  Related  Theory 

Since  UFS  is  heavily  based  on  kinetic  theory  and  the  Boltzmann  equation,  this  chapter 
briefly  introduces  the  kinetic  theory  of  gases,  the  Boltzmann  equation,  and  a  useful 
approximation  to  the  Boltzmann  equation.  Also,  the  numerical  schemes  for  the  solvers 
implemented  in  UFS  are  presented  here  to  give  the  reader  a  better  understanding  of 
how  UFS  performs  its  calculations. 

2.1  Kinetic  Theory  of  Gases 

Kinetic  theory  seeks  to  explain  macroscopic  flow  phenomena  by  investigating 
molecular  interactions.  For  example,  the  transport  phenomena  of  viscosity,  heat  con¬ 
duction,  and  diffusion  can  be  explained  by  the  molecular  transport  of  momentum, 
energy,  and  mass,  respectively. 

To  begin  with,  assume  a  molecule  can  be  modeled  as  a  rigid  sphere  with  no  in¬ 
ternal  structure.  The  only  physical  characteristics  of  the  molecule  are  then  an  effective 
diameter  and  a  mass.  This  molecular  model  is,  of  course,  a  simplification.  In  reality, 
electron  clouds  are  not  always  spherically  symmetric,  and  polyatomic  molecules  may 
have  an  intricate  structure  with  rotational,  vibrational,  and  electronic  energies.  How¬ 
ever,  the  rigid  sphere  model  is  a  good  starting  point.  Normally,  molecules  interact 
not  just  during  physical  “contact”  because  they  have  an  intermolecular  force  that  is 
a  function  of  their  distance  from  each  other.  Strangely,  this  force  is  repulsive  at  short 
distances,  but  attractive  at  large  distances,  while  the  more  familiar  gravitational  force 
is  solely  attractive.  The  rigid  sphere  model  instead  assumes  that  at  some  effective 
diameter,  there  is  an  infinite  repulsive  force,  and  everywhere  else  the  intermolecular 
force  is  zero.  The  dynamics  of  rigid  sphere  molecules,  therefore,  act  as  moving  billiard 
balls  on  a  table,  knocking  each  other  around  only  as  they  collide.  The  intermolecu¬ 
lar  force  can  be  modified  beyond  this  simplification  within  kinetic  theory,  as  will  be 
discussed. 

Each  of  the  molecules  has  a  translational  (or  kinetic)  energy,  which  is  dependent 
on  the  macroscopic  temperature.  However,  every  molecule  does  not  travel  in  space 
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with  the  same  speed,  nor  in  the  same  direction.  Some  molecules  have  a  higher  velocity 
than  the  mean  velocity,  while  some  have  a  lower  velocity  than  the  average.  Over  time, 
a  certain  molecule’s  velocity  will  change  due  to  collisions,  but  the  average  velocity  of  all 
the  molecules  per  unit  volume  will  remain  the  same  when  the  system  is  in  equilibrium. 
Because  there  are  so  many  molecules  per  unit  volume,  the  distribution  of  velocities 
over  a  unit  volume  in  velocity  space  can  be  represented  by  the  VDF.  Velocity  space 
has  as  its  axes  the  velocities  in  the  coordinate  directions  (Ci,  C 2 ,  and  C'3 ) ,  instead 
of  the  spatial  axes  (x,  y,  and  z).  Therefore,  regardless  of  where  a  molecule  lies  in 
physical  space,  its  point  in  velocity  space  corresponds  to  its  u-,  v-,  and  re-components 
of  velocity.  The  VDF  is  a  normalized  distribution  function,  meaning  that  the  volume 
under  its  curve  is  equal  to  unity: 

/OO 

f(Ci)dCi  =  1  (2) 

-00 

where  /  is  the  normalized  VDF,  C,  is  the  molecular  velocity  vector  (m/s),  and  dCi  = 
dC\  dCzdCs  is  an  element  in  velocity  space  (m3).  Also,  the  average  of  any  quantity 
Q  that  is  also  a  function  of  velocity,  can  be  calculated  from 

/OO 

QiCJfiCJdCi  (3) 

-00 

For  example,  the  average  velocity  in  the  x-direction  C\  is 

/OO 

C\f{Ci)dCi  (4) 

-00 

Equations  (2)  and  (3)  are  the  zeroth  and  first  moments  of  the  VDF,  respectively, 
where  the  nth  moment  of  the  VDF  Mn  is  defined  as 

/OO 

Q(Ci)f(Ci)dCi  (5) 
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The  VDF  for  a  gas  at  equilibrium  has  been  well  established,  and  was  first 


published  by  Maxwell  in  1860.  It  bears  his  name  as  the  Maxwellian  distribution: 


where  m  is  the  molecular  mass  (kg),  k  is  the  Boltzmann  constant  (1.38  x  10-23  kg  • 
m2/s2  •  K),  and  T  is  the  temperature  (K).  Following  the  development  by  Vincenti  and 
Kruger,  define  a  new  function  $  [37] : 


and  then 


feq(Q)  =  $(C1)$(C'2)$(C'3)  (8) 

In  other  words,  each  of  the  velocity  component  distribution  functions  $((7*)  is  statis¬ 
tically  independent  from  the  others.  The  Maxwellian  distribution  is  shown  in  Figure  2 
with  an  average  molecular  velocity  of  zero.  For  equilibrium,  the  distribution  is  sym¬ 
metric  about  its  mean.  Also,  for  all  VDFs,  the  macroscopic  velocity  corresponds  to 
the  average  velocity.  The  above  analysis  uses  the  thermal  velocity,  Ct  (m/s): 

Ci  —  Ci  —  ^  (9) 

where  q  is  the  molecular  velocity  (m/s)  as  seen  from  a  stationary  observer,  and 
Ci  is  the  average  molecular  (or  macroscopic  gas)  velocity  (m/s)  also  observed  by  a 
stationary  observer.  feq  can  then  be  written  as 


since  C2  =  (72  +  Cf  +  Cf  =  (|c*  —  c*|)2  .  So,  Figure  2  is  actually  showing  the 

thermal  velocity  distribution.  To  include  the  macroscopic  velocity,  the  distribution 
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Figure  2:  Non-dimensional  Maxwellian  distribution  function  for  only  one  component 
of  the  thermal  velocity  (centered  on  the  average  flow  velocity). 

in  Figure  2  needs  only  to  be  shifted  on  the  horizontal  axis.  As  noted  by  Vincenti 
and  Kruger,  in  a  dimensional  plot  of  $  versus  C'i,  increasing  only  the  temperature 
widens  the  curve  and  lowers  the  maximum,  while  increasing  only  the  molecular  mass 
has  just  the  opposite  effect  [37].  Since  the  full  VDF  is  actually  in  three-dimensional 
velocity  space,  it  is  difficult  to  visualize.  However,  Figure  3  gives  a  two-dimensional 
representation  of  the  equilibrium  VDF. 

The  VDF  is  also  used  to  calculate  other  macroscopic  flow  variables  such  as  p 
(kg/m3),  velocity  q,  temperature  T,  pressure  p  (N/m2),  shear  stress  TtJ  (N/m2),  and 
heat  flux  q,  (W/m2),  using  its  moments  [36]: 
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Figure  3:  Non-dimensional  Maxwellian  distribution  function  for  two  components  of 
the  thermal  velocity  (centered  on  the  average  flow  velocity). 


poo 

p  =  m  nf(ci)dci 

J  —  OO 

(ii) 

poo 

Ci  =  J  Cif(ci)dCi 

(12) 

1 

(13) 

1  f°° 

P=^P  (q  -  Ci)2f(ci)dCi 

^  J — OO 

(14) 

/»oo 

Tij  =  p  I  (Ci  -  Ci)(Cj  -  Cj)f(Ci)dCi 

(15) 

1  f°° 

Qi=  ~P  ( ^  -  Ci){cj  -  Cj)2 f  (ci)dci 

^  J — OO 

(16) 

where  R  is  the  specific  gas  constant  (J /kg-K),  and  n  is  the  number  density  (or  number 
of  molecules  per  unit  volume,  1/m3).1 


1It  is  interesting  to  note  that  temperature  and  pressure  are  directly  proportional  to  the  variance 

of  the  VDF  f  (a  -  Ci)2f(ci)dci. 
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Chapman  and  Enskog  independently  analyzed  small  departures  from  nonequi¬ 
librium  using  an  expansion  of  /  about  feq: 

f  —  feq(  1  +  $1  +  $2  +  $3  +  $4  +  •  •  • )  (17) 

where  d>  =  $(77n)  ,  and  <!>„  denote  higher-order  terms.  It  can  be  shown  that  when 
only  the  zeroth-order  term  is  considered  (the  first  term  in  the  expansion),  the  Euler 
equations  are  realized,  and  when  first-order  terms  are  included,  the  NS  equations  are 
achieved.  Higher-order  terms  have  been  included  to  develop  the  Burnett  (second- 
order  terms),  the  Super-Burnett  (third-order  terms),  and  the  Super-Super-Burnett 
(fourth-order  terms)  equations,  even  though  their  usefulness  seems  to  be  limited.2 
The  Chapman-Enskog  expansion  provides  a  useful  interpretation  of  the  Euler  and  NS 
equations.  The  Euler  equations,  because  they  only  include  the  zeroth-order  term, 
assume  the  flow  is  always  at  equilibrium,  showing  that  viscosity  and  heat  fluxes  are 
nonequilibrium  phenomena.  The  NS  equations,  because  they  include  at  most  the  first- 
order  terms,  assume  the  flow  only  experiences  small  departures  from  nonequilibrium. 
Therefore,  continuum  relations  become  invalid  for  flows  that  depart  significantly  from 
nonequilibrium. 

Since  it  is  very  expensive  computationally  to  follow  every  particle  in  a  flow,  it 
would  be  useful  to  solve  for  the  VDF  directly  using  some  governing  equation.  The 
Boltzmann  equation  is  such  a  governing  equation  for  the  VDF,  describing  how  the 
VDF  evolves  in  space  and  time.  The  number  of  molecules  of  class  c,  (meaning  with 
velocity  q)  in  a  physical  space  volume  element  dVx=  dx  dy  dz  (m3)  and  in  a  velocity 
space  volume  element  dVc=  dcidczdc?,  (m3/s3)  is  nf(ci)dVxdVc,  where  n  is  the 
number  of  molecules  per  unit  physical  volume  (1/m3).  The  number  of  molecules  in 

2For  a  study  on  one-dimensional  shock  structure  in  a  monatomic  gas  using  the  Super-Burnett 
and  Super-Super-Burnett  equations,  see  [35]. 
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these  volume  elements,  over  time,  is  then 


^[nf(ci)]dVxdVc  (18) 

which  can  only  be  a  result  of  convection  of  molecules  in  physical  space,  convection 
of  molecules  in  velocity  space,  or  collisions  within  dVx  which  deplete  or  increase  the 
number  of  molecules  of  class  q.  Considering  all  six  sides  of  the  volume  element  dVx, 
the  convection  in  physical  space  is  modeled  as 

-cj-^lnfic^dV.dV.  (19) 

The  convection  in  velocity  space  would  be  caused  by  an  acceleration  F*  (m/s2)  due 
to  some  external  force,  such  as  gravity  or  an  electromagnetic  field: 

-  -^[Fjnf(ci))dVxdVc  (20) 

The  general  term 

{J^/te)]}  dvxdvc  (21) 

is  used  to  describe  how  molecules  of  class  c*  increase  due  to  collisions.  Combining 
these  terms  and  dividing  by  dVx  dVc ,  the  full  Boltzmann  equation  is 

JWte)]  +  cj-Srinf(ci)}  +  4r[Finf(ci)\  =  { JWte)]  j  (22) 

The  term  on  the  right-hand  side  is  called  the  collision  term,  or  the  collision  inte¬ 
gral,  because  its  direct  calculation  involves  multi-dimensional  integration.  Vincenti 
and  Kruger  develop  the  collision  integral  for  hard-sphere  molecules  (with  a  spherical 
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intermolecular  potential),  given  as  [37] 

C  n  ^  r oo  p2tt  /*7t/2 

\  \  =  /  /  /  [/(4)/(C0  -  /(q)/(C*)]^2  sin^  COS  xpdxpdedV^ 

lOZ  J  coll  J-ooJO  j 0 

(23) 

where  q  and  Q  are  the  pre-collisional  velocities  and  c'  and  ('  are  the  post-collisional 
velocities  of  the  two  molecules  (m/s),  g=  |</t  —  q|  is  the  relative  molecular  speed 
(m/s),  d  is  the  diameter  of  the  molecules  (m),  d2  sin  ip  cos  ipdipde  is  the  collisional 
differential  cross-section  on  a  sphere,  and  dV(  is  the  volume  element  <7(/  d/2  d(3  in 
velocity  space.  The  first  term  in  the  collision  integral  represents  inverse  collisions 
which  replenish  the  number  of  molecules  of  class  q ,  while  the  second  term  represents 
collisions  which  deplete  the  number  of  molecules  of  class  q.3 

A  couple  of  observations  can  be  made  at  this  point  about  the  Boltzmann  equa¬ 
tion.  The  first  observation  is  that  the  Boltzmann  equation  is  an  integro-differential 
equation,  meaning  that  the  VDF  is  the  argument  of  both  integrals  and  differentials, 
making  it  very  difficult  to  evaluate.  The  second  is  that  the  Boltzmann  equation 
assumes  that  only  bimolecular  collisions,  or  collisions  involving  only  two  molecules, 
occur.  Collisions  involving  more  than  two  molecules  are  more  likely  to  occur  when 
the  density  is  high.  However,  the  densities  encountered  in  this  work  are  low  enough 
that  this  assumption  applies. 

The  Boltzmann  equation  can  be  arrived  at  from  the  more  general  Liouville 
equation,  which  is  in  6A^-dimensions  describing  how  every  possible  microstate  (defined 
by  the  position  and  momentum  of  every  particle)  is  tracked  in  phase  space  (the  space  of 
all  states  at  which  the  system  can  exist).4  The  Boltzmann  equation  can  be  rigorously 
derived  from  the  Liouville  equation  by  only  considering  the  one-particle  distribution 
function  /  and  integrating  over  all  dimensions  except  for  q  and  :q.  An  in-depth 

3 For  identical  incidence  angles,  the  velocities  after  a  deplenishing  collision  are  equal  to  c,  and  Q , 
while  the  the  velocities  after  a  replenishing  inverse  collision  are  equal  to  ct  and  Q. 

4Here,  N  denotes  the  number  of  molecules. 
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discussion  of  this  subject  is  beyond  the  scope  of  this  work,  and  the  reader  is  referred 
to  the  book  by  Harris  for  a  more  complete  introduction  [14]. 

Even  though  the  Boltzmann  equation  has  been  around  for  many  years,  an  an¬ 
alytic  solution  has  not  been  found.  In  an  attempt  to  avoid  the  costly  calculation  of 
the  collision  integral,  approximate  models  to  the  Boltzmann  equation  have  been  pro¬ 
posed.  One  such  model  was  developed  by  Bhatnagar,  Gross,  and  Krook,  called  the 
Bhatnagar-Gross-Krook  (BGK)  collision  model  [37].  The  BGK  model,  even  though 
it  involves  a  simplification  of  the  collision  term,  tries  to  retain  some  features  of  the 
same.  The  collision  term  is  replaced  by 

{|[»/ta)]}  =Mr-f)  (24) 

where  u  is  a  collision  frequency  (Hz,  or  1/s),  feq  is  the  equilibrium  (or  Maxwellian) 
distribution,  and  /  is  the  variable  VDF.  The  full  BGK  model  is  therefore 

+  G^-[n/(o)]  +  [*>/(<*)]  =  np(yfeq  ~  /)  (25) 

Since  v  is  a  collision  frequency,  it  could  also  be  written  as  1/iy,  where  rv  is  a  local 
relaxation  time  (s).  Therefore,  the  BGK  collision  term  takes  on  the  same  form  as  that 
used  in  the  vibrational  relaxation  equation.  In  light  of  this  relationship,  the  BGK 
model  essentially  describes  how  a  gas  in  nonequilibrium  relaxes  (or  equilibrates)  to 
a  state  of  equilibrium.  So,  one  of  the  drawbacks  of  the  BGK  model  is  that  it  loses 
accuracy  or  validity  as  the  departure  from  nonequilibrium  increases.5  Vincenti  and 
Kruger  also  note  that  since  feq  is  a  function  of  c,  and  T,  and  since  c,  and  T  are 
calculated  as  moments  of  the  VDF  (see  Equations  (12)  and  (13)),  the  BGK  equation 
is  still  a  nonlinear  integro-differential  equation  [37]  ,6 

5Xu  and  Guo  extended  the  BGK  model  further  into  the  nonequilibrium  regime  by  replacing  the 
one-stage  BGK  collision  model  with  a  two-stage  model  [40]. 

6See  Equation  (10). 
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Another  aspect  of  the  molecular  model  that  should  be  considered  is  the  inter- 
molecular  potential.  Molecules  far  apart  tend  to  weakly  attract  one  another,  whereas 
molecules  close  in  proximity  tend  to  strongly  repel  one  another.  The  Hard  Sphere 
(HS)  model,  because  it  does  not  account  for  any  other  interaction  besides  collisions, 
assumes  only  an  infinitely  large  repulsive  force  during  a  collision.  The  Sutherland 
model,  on  the  other  hand,  also  allows  for  an  attractive  force  that  decreases  with 
distance.  When  solving  the  Boltzmann  equation,  one  usually  must  define  which  in- 
termolecular  potential  model  one  is  using.  With  the  BGK  model,  however,  since  the 
collision  integral  is  replaced  entirely,  no  explicit  potential  model  is  used.  Rather,  a 
certain  collision  model  is  specified  in  the  calculation  of  v  (this  study  uses  the  HS 
collision  model). 

2.2  UFS  Methods 

UFS  is  built  on  the  Gerris  code  [2],  Gerris  is  an  open-source  code  which  was 
initially  released  in  2005  and  is  supported  by  New  Zealand’s  NIWA  (National  Institute 
of  Water  and  Atmospheric  research).  It  is  a  time-dependent,  second-order  (in  time  and 
space),  finite- volume  solver  for  incompressible  flows,  with  the  option  of  using  either  an 
Euler  or  a  NS  scheme,  and  is  capable  of  parallel  computing  using  the  MPI  (Message- 
Passing  Interface)  library.  It  uses  Cartesian  quadtree  (two-dimensional)  or  octree 
(three-dimensional)  grids,  meaning  that  all  cells  are  either  squares  or  cubes,  whose 
dimensions  only  differ  by  some  factor  of  two.  The  grid  structure  will  be  discussed 
more  in  Chapter  III  in  reference  to  the  particular  grids  used  in  this  study. 

UFS  utilizes  both  continuum  and  kinetic  solvers  when  computing  flow  parame¬ 
ters.  The  user  can  specify  which  solvers  to  use,  and  where  in  the  domain  they  should 
be  called.  The  kinetic  solvers  are  either  approximations  of  the  Boltzmann  equation 
(such  as  the  BGK  model)  or  DNS  (Direct  Numerical  Simulation)  of  the  collision  in¬ 
tegral  using  some  intermolecular  force  model.  As  of  2006,  four  intermolecular  force 
models  have  been  included  in  UFS,  namely  the  HS  model,  the  inverse  power  repulsive 
potential,  the  Lennard- Jones  potential,  and  the  Coulomb  potential  [21]. 
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As  explained  by  Kolobov  et  al.,  the  Boltzmann  equation  as  Equation  (22)  can 
also  be  framed  in  the  following  form  when  neglecting  external  forces  [21]: 

%  +  V,  •  (if)  =  Hf,  f)  (26) 

where  r  is  a  position  vector  in  physical  space  (m),  £  is  the  velocity  vector  (m/s), 
and  /(/,  /)  is  the  collision  integral.  In  order  to  solve  Equation  (26)  numerically,  a 
Cartesian  mesh  is  created  in  velocity  space  with  nodes  £$  and  cell  size  (m/s).7 
It  is  claimed  that  by  using  this  mesh,  Equation  (26)  is  reduced  to  a  system  of  linear 
equations  (hyperbolic)  in  physical  space  with  the  nonlinear  source  term 

=  (27) 

where  the  subscript  i  indicates  that  the  parameter  is  evaluated  at  nodes  i.  Since  UFS 
is  built  on  the  Gerris  solver,  which  utilizes  only  Cartesian  grids  in  physical  space, 
the  computational  grid  in  physical  space  in  UFS  is  also  Cartesian.  Equation  (27) 
is  solved  in  two  stages — collisionless  flow  and  relaxation,  divided  by  an  intermediate 
time  stage. 

For  the  collisionless  flow  stage,  the  collision  integral  is  assumed  to  be  zero,  and 
an  explicit  finite  volume  scheme  is  used: 

£*k  _  rk—1 

V  At  f  ‘  n)face  -4 face  ^face  =  0  (28) 

face 

where  j  is  the  cell  number  in  physical  space,  k  is  the  time  index,  *  indicates  an 

intermediate  time  level,  n  is  the  unit  outward  normal  vector  to  the  cell  face,  V  is  the 

"71  is  the  value  of  the 
face 

VDF  on  the  cell  face.  To  calculate  fk^ce ,  interpolation  schemes  are  used  (either  first- 
or  second-order),  with  three  options  for  the  second-order  limiter:  none,  minmod,  or 

7In  future  work,  CFDRC  hopes  to  include  automatic  mesh  refinement  in  velocity  space,  which 
would  save  computational  expense  as  well  as  make  velocity-space  grid  independence  studies  simpler. 


cell  volume  (nr),  Sface  is  the  face  surface  area  (nr),  and  / 
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van  Leer.  The  minmod  limiter,  which  is  the  most  conservative  of  the  second-order 
limiters,  is  used  in  this  study. 

For  the  relaxation  stage, 


fk  _  f*k 
•J  ij  J  ij 

At 


-u*k  f*k  +  $*.' 

IJ  J  IJ  1  11 


IJ 


(29) 


where  v  is  the  collision  frequency  (how  often  molecules  in  class  £  are  depleted)  and 
<F  is  the  inverse  collision  integral.  Similar  to  the  BGK  model,  UFS  replaces  here  the 
complex  collision  integral  with  two  terms  {—vf  and  <F),  illustrating  again  that  the 
BGK  model  describes  how  the  distribution  function  relaxes  to  equilibrium. 

The  continuum  solvers  are  also  based  on  the  Boltzmann  equation  because  macro¬ 
scopic  flow  parameters  can  be  obtained  through  the  moments  of  the  VDF.  Also,  cou¬ 
pling  of  the  continuum  and  kinetic  solvers  is  easier  when  they  are  all  based  on  the 
Boltzmann  equation.  Again,  Kolobov  et  al.  present  the  following  material  concerning 
the  Euler  and  NS  solvers,  which  is  included  here  for  completeness  [21]. 

The  Euler  equations  can  be  presented  in  the  form 


dY  OF  dG  dH_ 

dt  +  dx  dy  +  dz 

where 

Y  =  [ p ,  pu,  pv,  pw ,  E] 

F  =  [pu,  p/2  +  pu2,  puv,  puw,  u{E  +  p)] 
G  =  [pv,  puv,  p/2  +  pv2,  pvw,  v(E  +  p)] 
H  =  [pw,  puw,  pvw,  p/2  +  pur ,  w(E  +  /.j] 


(30) 


(31) 

(32) 

(33) 

(34) 
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and  p  is  the  density  (defined  as  p  =  mn  ),  u,  v,  and  w  are  the  flow  velocities  in  the 
x-,  y- ,  and  ^-directions  (m/s),  respectively,  E  is  the  energy  (J): 


E 


3 

2 


pT  +  p(u 2  +  v2  +  w2) 


(35) 


T  is  the  temperature,  and  p  =  nkT  is  the  pressure.  Discretizing  Equation  (30)  with 
an  explicit  finite  volume  scheme,  and  moving  spatial  derivatives  to  the  right-hand 
side, 


yrn+1 
1  ijk 


_  -yrn 
1  ijk 


At 

i+l/2,j,k 


jpn 


T7in 

^  i—l/2,j,k 


Ax 


+ 


fin 

L*i,j+l/2  ,k 


fin 

^ i, j-l/2, k 


Ay 


+ 


TTn 

a  i,j,k+ 1/2 


-rrn 

ij,k— 1/2 


Az 


(36) 


where  n  denotes  the  time  step,  i,  j,  and  k  denote  the  cell  nodes  in  physical  space 
for  the  x-,  y-,  and  ^-directions,  Yr-jk  is  the  cell-averaged  value  at  time  step  n,  and 
F?+1/2J>k,  and  HE^k+l/ 2  indicate  the  fluxes  on  the  cell  faces  along  the  x-, 

y-,  and  ^-directions,  respectively.  The  fluxes  are  calculated  using  moments  of  the 
distribution  function: 


1  ftn+1  f 

Fi±i/2j,k  = -r;  /  il>€xf{xi±i/2,yj,zk,t,g)d£dt  (37) 

Jtn  Jl3 

1  ftn+1  f 

Gij±i/2,k  =  ~r;  /  /  ‘•l>€yf{xi,yj±i/2,Zk,t,g)d£dt  (38) 

Jtn  Jr3 

1  ftn+1  f 

HiJtk±1/2  =  —  /  yjy  zk±1/2,  t,  £)d£dt  (39) 

Jtn  Jr3 


where  M3  indicates  that  the  integral  is  over  all  real  values  in  velocity  space,  and  -0 
signifies  the  collisional  invariants,  given  as  -?/>(£)  =  [1,  £2]  .  The  collisional 
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invariants  are  so  named  because  they  have  the  property  that8 


[  il>I(f,f)d£  =  0 

Jr3 

The  VDF  at  the  cell  faces  is  calculated  as 


(40) 


f(x j+1/2,  Vi,  zt,  t,  Z)  =  H[(x]f l"  +  (1  -  //[{«])/, 


eq 

r 


(41) 


where  f[q  and  are  the  equilibrium  VDFs  on  the  cell’s  left  and  right  faces,  respec¬ 
tively: 


if 


Pi- 1/2 


(jTT"- 1/2 ) 


\  3/2 


exp 


fei  = _ ^+1/2 _ exD 

Jr  *  x  3/2 


nTn 

l 


+1/2  J 


(£x  1/2, j,k)  T  (Cy  1/2, j, k)  T  (£z  1/2, j,k) 

rpn 

1i-l/2 

(42) 

(&e  ~  lh+l/2,j,fc)'~  ~F  ~  ^+l/2,j,fc)  T  (£z  ~  ‘^h+l/2,j,fc) 

r Tn 

±i+ 1/2 

(43) 


and  i7[£]  is  the  Heaviside  step  function: 


|  1,  £>0 
( o,  e  <  o 


(44) 


Hence,  the  Euler  scheme  uses  only  Maxwellian  VDFs,  which  is  consistent  with  the 
zeroth-order  Chapman-Enskog  expansion.  If  a  Boltzmann  solver  is  being  used  in  a 
neighboring  cell,  then  the  parameters  required  for  the  calculation  of  in  Equa¬ 
tions  (42)  and  (43)  are  found  from  moments  of  the  VDF  in  the  neighboring  cell. 
Similarly,  a  Boltzmann  cell  would  assume  a  Maxwellian  VDF  in  a  neighboring  Euler 
cell.  The  first-order  Euler  scheme  uses  only  a  two-cell  stencil  {x^j  and  xi+ij,  or  Xi-ij 

8For  a  more  detailed  description  of  the  collisional  invariants,  the  reader  is  referred  to  Sone’s  book 
on  kinetic  theory  [36]. 
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and  Xij),  while  the  second-order  discretization  uses  a  three-cell  stencil  (aq_ ij,  xUj , 
and  Xi+ij)  with  one  of  the  three  available  limiters. 

The  kinetic  NS  solver  is  also  based  on  the  solution  of  the  Boltzmann  equation, 
although  in  a  somewhat  different  fashion  [39].  It  is  taken  from  the  BGK  model: 


df 

dt 


+  Vr-(£/) 


/eg-/ 

r 


(45) 


where  r  is  the  intercollision  relaxation  time  (s)  defined  as  r  =  y/p  ,  and  p  is  the 
dynamic  viscosity  (N-s/m2).  A  directional  splitting  method  is  used  to  reduce  these 
multi-dimensional  equations  to  a  set  of  one-dimensional  equations,  and  then  solved 
analytically  [39].  For  a  one-dimensional  case,  the  VDF  is  solved  to  be 

f(x,  £,  t )  =  e~t,T  f0(x  -  £,  0)  +  -  [  feq(xi,  £,  tt)  e~{t~tl)/r  dtt  (46) 

r  Jo 

where  ti  is  a  dummy  integration  variable  for  time  (s),  X[  is  the  trajectory  of  the 
particles  (m),  calculated  as  xi  =  x  —  u(t  —  ti)  ,  and  the  functional  dependencies  on 
y  and  z  have  been  omitted.  The  VDF  is  calculated  for  each  cell  i  with  faces  i  +1/2 
and  i  —  1/2,  with  the  VDF  at  some  initial  moment  fo : 


fo  =  fix,  £,t  =  0) 

=  //69[1  +  aix  ~  Tiaiix  +  ^-/)](1  —  H[x])  +  //9[  1  +  arx  —  r(ar^x  +  Ar)]H[x]  (47) 
and  the  equilibrium  VDF  feq  is  calculated  from 

feq  =  /0e9[  1  +  (1  -  H[x])atx  +  H[x]drx  +  At]  (48) 
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where 


al,r  =  al,r  +  °h,Xx  +  af,r£y  +  atr£z  +  at,r((x  +  £y  +  £z) 

^l,r  =  +  <5f,r6r  +  +  ^t,r^z  +  bf;r(£x  +  +  £z) 

Ar  =  Af>r  +  A  lr£x  +  Azlr^y  +  Af^z  +  Afr(£2  +  £y  +  Cf) 


A  =  A1  +  a2^  +  A%  +  a4^  +  A5(^  +  %  +  42) 


(49) 

(50) 

(51) 

(52) 


are  polynomial  functions  of  the  local  constants  a\  r,  a]r,  Allr,  and  A1  (  i  —  1,  2,  3,  4,  5  ), 
and  /g9  =  feq(x,  £,  t  =  0)  ,  which  are  all  determined  by  the  method  given  by 

Li  et  al.  [23] 

Nonequilibrium  is  accounted  for  in  the  VDF  with  the  terms  r//9(a;C-  +  Ai ) 
and  rf^q(ar^x  +  Ar )  in  Equation  (47).  As  long  as  xfiqr{ai^x  +  A^r)  <C  1  ,  the 
approximation  given  here  for  the  VDF  is  consistent  with  the  first-order  Chapman- 
Enskog  expansion,  and  therefore  recovers  the  NS  equations. 

If  a  NS  cell  has  a  neighbor  Boltzmann  cell,  then  the  NS  cell  is  given  a  velocity 
grid  identical  to  the  neighboring  cell,  and  the  VDF  f0  =  //9[  1  —  r(a/ir^n  +  Aitr)\  is 
created  on  the  cell  interface  where  £n  is  the  normal  velocity  to  the  cell  face  (m/s).  /// 
is  calculated  using  the  macroparameters  from  the  NS  cell,  and  the  coefficients  for  a^r 
and  ALr  are  obtained  using  gradients  of  macroparameters  from  both  cells. 

The  CFL  (Courant-Friedrichs-Lewy)  condition  is  used  to  determine  the  time 
step  for  both  the  Euler  and  NS  solvers.  Since  all  of  the  simulations  presented  in  this 
work  are  unsteady,  each  cell  uses  the  same  time  step  (the  minimum  time  step  for  the 
entire  domain).  The  CFL  condition  is 


A  t  = 


CFL  x  h 

max(|[/  +  3VT|,  \U  —  3\/T|) 


(53) 


where  At  is  the  non-dimensional  time  step,  h  is  the  local  non-dimensional  cell  size, 
U  is  the  non-dimensional  flow  velocity  defined  as  U2  =  u2  +  v2  +  w2  ,  T  is  the 
non-dimensional  temperature,  and  max(| U  +  3\/T|,  | U  —  3\ZT|)  =  |£max|  for  the 
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Boltzmann  solver  with  £ max  being  the  non-dimensional  molecular  velocity.  The  CFL 
is  typically  set  to  0.5.  The  time  required  to  run  each  simulation  is  reported  in  Chap¬ 
ter  IV. 

A  well-known  consequence  of  the  BGK  model  is  that  the  Prandtl  number  can 
only  be  unity  [37]: 


(54) 


where  cp  is  the  specific  heat  at  constant  pressure  (J/kg-K),  and  k  is  the  coefficient 
of  thermal  conductivity  (W/m-K).  The  Prandtl  number  is  the  ratio  of  the  viscous 
diffusion  rate  to  the  thermal  diffusion  rate,  and  provides  an  indication  of  the  rela¬ 
tive  importance  of  each.  Thermal  diffusion  considers  the  transfer  of  energy  through 
conduction,  while  viscous  diffusion  refers  to  the  transfer  of  energy  due  to  molecular 
mixing.  When  using  the  BGK  model  (including  the  kinetic-based  NS  solver),  one 
can  accurately  model  either  viscosity  /i  or  the  coefficient  of  thermal  conductivity  k, 
but  not  both.  Xu  presents  a  Prandtl  number  fix  based  on  altering  the  heat  flux  Q 
(W/m2)  (see  Equation  (16))  [39].  The  energy  flux  (the  last  term  in  the  vectors  F . 
G,  and  H)  is  modified  by  adding  the  correction  term  (Pr~l  —  1  )Q  to  it,  where  Q  is 
calculated  on  the  cell  faces  using  polynomial  interpolation  of  the  VDF,  and  Pr  is  a 
variable  Prandtl  number.  Results  are  reported  for  shockwave  structure,  showing  that 
the  regular  NS  solver  (based  on  fluid  continuum  assumptions)  and  the  kinetic-based 
NS  solver  show  good  agreement  [24,39]. 

Since  UFS  is  a  hybrid  code,  it  has  kinetic  solvers  in  addition  to  the  kinetic- 
based  continuum  solvers.  The  kinetic  solver  used  in  this  work  is  the  BGK  model. 
However,  since  the  BGK  model  does  not  calculate  a  collision  integral,  it  does  not 
inherently  have  the  capability  to  account  for  internal  molecular  energies,  and  is  only 
used  for  the  Argon  simulations.  In  order  to  retain  the  computational  efficiency  of  the 
BGK  model,  UFS  uses  a  three-temperature  BGK  (3T-BGK)  model  for  flows  involving 
internal  energies  (utilized  for  the  V2  simulations).  Most  of  the  details  of  this  model 
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are  not  available  for  public  release  since  UFS  is  regulated  by  ITAR  (International 
Traffic  in  Arms  Regulations).  However,  as  a  brief  introduction,  the  3T-BGK  model 
assumes  that  there  are  three  different  temperatures  Ttr  (K)  (translational),  Trot  (K) 
(rotational),  and  Tvib  (K)  (vibrational).9  The  equilibrium  temperature  Teq  (K)  can 
be  obtained  from  a  weighted  average  of  the  temperatures: 

rpeq  _  +  KrTrot  +  KvTvib 

3  +  Kr  +  K^q 

where  Kr  and  Kv  are  the  rotational  and  vibrational  degrees  of  freedom,  respectively, 
and  K^q  =  Kv(Teq )  ,10  The  Maxwellian  VDF  feq  is  then  calculated  as 

r  =  feq(Ttr)feq(Trot)rq(Tvib )  (56) 

which  assumes  that  the  VDFs  for  each  of  the  three  temperatures  are  statistically 
independent  of  one  another. 

In  order  to  successfully  couple  a  continuum  solver  with  a  kinetic  solver,  a  switch¬ 
ing  parameter  is  required.  Such  a  parameter  could  indicate  if  the  continuum  as¬ 
sumption  is  appropriate  for  any  given  cell.  Hence,  the  switching  parameter  is  also 
sometimes  referred  to  as  the  continuum  breakdown  parameter.  As  mentioned  by 
Kolobov  et  al.,  it  is  important  to  use  an  appropriate  switching  parameter  in  order  to 
avoid  negative  values  in  the  VDF  (which  is  non-physical)  [21].  For  example,  Schrock 
and  Carr  both  investigated  entropy  generation  as  a  means  of  quantifying  continuum 

9Josyula  et  al.  used  a  three-temperature  BGK  model  as  well,  but  the  three  temperatures  were 
instead  translational  temperatures  Tx,  T,n  and  Tz  [19].  Also,  Josyula  et  al.  implemented  a  two- 
temperature  BGK  model  for  the  translational  and  rotational  temperatures,  Tb  and  Trot  [18]. 

10  N2  has  two  (not  three)  rotational  degrees  of  freedom  (  Kr  =  2  ).  If  the  N2  molecule  is  modeled 
as  a  dumbbell  (or  rigid  rotor),  then  the  energy  of  rotation  about  the  axis  which  passes  through  the 
connector  is  very  small  compared  to  the  energy  of  rotation  through  the  other  two  principal  axes.  N2 
has  three  vibrational  degrees  of  freedom  (  Kv  =  3  ). 


(55) 
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breakdown  [7, 34] .  The  switching  parameter  used  in  the  first  part  of  this  work  is 


where  all  values  are  non-dimensional,  and  U2  =  u2  +  v2  +  w2  ,  has  been  shown  by 
Kolobov  et  al.  to  correctly  couple  the  solvers  near  a  shockwave  at  moderate  Knud- 
sen  numbers  [21].  A  user-specified  threshold  value  for  Sns  is  provided  in  an  input 
file,  and  Sns  is  calculated  for  each  cell.  If  a  cell’s  value  for  Sns  is  larger  than  the 
threshold,  then  that  cell  is  flagged  as  a  Boltzmann  cell.  Otherwise,  it  is  flagged  as 
a  continuum  cell.  Thus,  by  decreasing  the  Sns  threshold  value,  one  decreases  the 
number  of  continuum  cells  in  the  domain,  while  increasing  the  threshold  value  makes 
the  Boltzmann  region(s)  smaller. 

The  threshold  value  used  for  the  entropy-shock  interaction  study  is  based  on 
the  density  gradient  and  local  Knudsen  number  [21,32]: 

Sp  =  kJ—  (58) 

P 

since  the  entropy  spot  has  an  accompanying  density  fluctuation  profile. 
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III.  Methodology 


3. 1  Shock  Structure 

In  order  to  study  shock  structure  using  UFS,  it  was  determined  that  a  stationary 
two-dimensional  shock  should  be  used.  It  was  thought  that  perhaps  a  shocktube 
domain  should  be  used,  since  the  Alsmeyer  data  was  collected  in  a  shocktube  [3]. 
However,  the  increased  computational  cost  of  resolving  the  expansion  wave  and  the 
contact  surface,  as  well  as  the  expanded  domain  that  would  be  required  made  it 
undesirable.  Therefore,  the  shock  reference  frame  was  chosen  over  the  laboratory 
reference  frame.  Besides,  the  structure  of  the  shock  is  more  easily  observed  in  the 
shock  reference  frame,  which  remains  stationary.  Since  UFS  is  built  on  the  Gerris 
software,  which  uses  quadtree  hnite- volume  discretization  in  two-dimensional  physical 
space,  the  only  way  to  have  a  one-dimensional  simulation  is  to  require  that  all  of  the 
cells  have  the  same  dimensions.  If  one  is  only  interested  in  one-dimensional  shock 
structure,  then  this  grid  is  perhaps  a  good  option.  But,  a  goal  of  this  research  is 
to  first  study  shock  structure,  and  then  to  study  how  imperfections  in  the  upstream 
flow  (such  as  a  spot  of  entropy)  affect  both  the  shock  and  the  downstream  flow 
characteristics.  Such  a  study  requires  a  two-dimensional  physical  domain. 

Multiple  grids  (coarse,  medium,  and  fine)  were  used  to  determine  grid  conver¬ 
gence  in  physical  space.  Figure  4  shows  the  coarse,  medium,  and  fine  grids  for  the 
Argon  and  Nitrogen  simulations.  The  dimensions  of  the  large  cells  for  the  medium 
grids  are  A0  =  Ax0  =  A y0  =  Ai  (the  upstream  mean  free  path),  and  the  small 
cells  are  of  dimension  A2  =  A.r2  =  Ay2  =  Ai/4  ,  where  Ao  is  the  dimension  for 
base  cells  (m)  (refinement  Level  0)  and  A2  is  the  dimension  for  cells  of  refinement 
Level  2  (m).1  For  the  coarse  grids,  A0  =  2Ai  and  A2  =  Ai/2  .  For  the  fine 
grids,  A0  =  Ai/2  and  A2  =  A1/8  .  The  reference  length  Lref  (m),  or  the  length 
by  which  all  other  lengths  are  normalized,  is  Lref  =  Ao  for  each  of  the  grids.  The 

1UFS  is  built  on  the  Gerris  software,  which  uses  Cartesian  quadtree  (or  octree  for  three  di¬ 
mensions)  finite  volume  discretization.  All  cells  are  squares,  and  the  base  (or  largest)  cells  are 
considered  Level  0  (  Ao/ Lref  =  2~°  =  1  ).  More  refined  cells  have  higher  levels  of  refinement, 
with  An/Lref  =  2~n  . 
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(a)  Argon  (Coarse) 


(d)  Nitrogen  (Coarse) 


(f)  Nitrogen  (Fine) 

Figure  4:  Coarse,  medium,  and  fine  grids  in  physical  space. 

medium  grid  for  Argon  is  8OA1  x  Ai  with  the  standing  shock  located  halfway  between 
the  upstream  and  downstream  boundaries  (40Ai  away  from  both  boundaries),  where 
the  flow  is  from  left  to  right.  The  refined  area  is  14Ai  upstream  and  downstream  of 
the  shock,  respectively.  The  medium  grid  for  N2  is  8OA1  x  Ai.  The  shock  is  located 
20Ai  from  the  upstream  boundary,  with  refined  regions  7.5Ai  upstream  and  4.5Ai 
downstream  of  the  shock,  respectively.  In  all  cases,  the  medium  grids  were  found  to 
provide  sufficient  accuracy  (the  density  profiles  of  the  medium  and  fine  grids  were 
within  1%  at  every  point  for  both  the  Argon  and  Nitrogen  simulations),  and  were 
therefore  used  for  the  results  presented  in  Chapter  IV. 

A  value  of  Sns  =  0.001  was  chosen  for  the  Argon  and  Nitrogen  simulations, 
with  the  subsequent  domain  decomposition  shown  in  Figure  5  (cells  in  light  blue  are 
flagged  as  continuum  cells,  while  cells  in  yellow  are  flagged  as  Boltzmann  cells).  The 
Boltzmann  domain  upstream  of  the  shock  increases  with  increasing  Mach  number  due 
to  the  1/U2  term  in  Equation  (57),  while  the  Boltzmann  domain  downstream  of  the 
shock  decreases  due  to  the  predicted  thinning  of  the  shock  by  the  BGK  model.  This 
result  will  be  further  discussed  in  Chapter  IV.  The  domain  decomposition  for  the 
Nitrogen  cases  are  shown  in  Figure  6  for  the  same  threshold  value  Sns  =  0.001  and 
color  coding  as  in  Figure  5.  The  same  observations  can  also  be  made  here,  namely 


(k)  Mi  =  9.0 

Figure  5:  Visualizations  of  Boltzmann  (yellow)  and  continuum  (light  blue)  cells  in 
Argon  at  steady-state  with  Sns  =  0.001  .  The  shock  is  located  in  the 
middle  of  the  refined  grid  region. 

that  the  upstream  Boltzmann  domain  increases  with  increasing  Mach  number,  while 
the  downstream  Boltzmann  domain  decreases  due  to  the  predicted  thinning  of  the 
shock. 

A  convergence  study  was  performed  on  the  switching  parameter  Sns  to  deter¬ 
mine  if  Sns  —  0.001  accurately  reflects  the  full  Boltzmann  solution.  All  of  the  cases 
were  simulated  with  only  the  Boltzmann  solver  turned  on  (either  BGK  or  3T-BGK), 
with  no  more  than  a  1%  difference  in  density  values  between  them  and  the  results 
from  the  coupled  solvers. 
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(a)  Mi  =  1.53 
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(b)  Mi  =  1.7 
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(c)  Mi  =  2.0 


l  l  l  l  l  l  l  l  l  l  l l l l  l  l  l  l l l  l l l  l l l  l l  l ll  l  l l  l l l  l l  l  l  l  l  l  l  l  l  l  l  l  l  l  l  l  l 


(d)  Mi  =  2.4 
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(e)  Mi  =  2.8 
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(f)  Mi  =  3.8 


Figure  6:  Visualizations  of  Boltzmann  (yellow)  and  continuum  (light  blue)  cells  in 

at  steady-state  with  S^s  =  0.001  .  The  shock  is  located  in  the  middle  of 
the  rehned  grid  region. 


Also,  a  grid  independence  study  was  performed  in  velocity  space.  For  the  Argon 
simulations  on  the  medium  velocity  grid,  the  range  in  velocity  space  was  —8 Vref  to 
121 with  a  node  spacing  of  0.2 Vref  for  Mi  <  8  ,  and  —15 Vref  to  15Vref  with 
a  node  spacing  of  0.3 Vref  for  Mi  >  8  in  both  the  u-  and  u-directions  (100  x  100 
nodes),  where  Vref  is  the  reference  velocity  (m/s),  chosen  here  to  be  the  upstream 
non-dimensional  thermal  velocity,  defined  as  \/2RT\ ,  where  the  upstream  temperature 
7i  is  non-dimensionalized  by  the  reference  temperature  Tref  (K).  The  fine  and  coarse 
grids  had  twice  and  half  as  many  nodes  as  the  medium  grid,  respectively.  The  density 
values  differed  by  less  than  1%.  For  the  Nitrogen  simulations  on  the  medium  velocity 
grid,  the  range  in  velocity  space  varied  more  with  Mi  (±5 Vref  with  20  x  20  nodes 
for  Mi  <  2  ,  ±6 Vref  with  24  x  24  nodes  for  M\  =  2.4  ,  ±7 Vref  with  28  x  28  nodes 
for  Mi  =  2.8  ,  and  ±9V^.e/  with  36  x  36  nodes  for  Mi  =  3.8  ),  while  the  spacing 
was  kept  at  0.5 Vref  in  the  u-  and  u-directions.  The  fine  and  coarse  velocity  grids  had 
half  and  double  the  spacing  in  the  medium  grid,  respectively.  Density  values  for  the 
medium  grids  all  agreed  with  those  for  the  fine  grids  within  1%,  showing  sufficient 
convergence. 
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For  the  Argon  cases,  Mach  numbers  Mi  =  1.55  ,  1.75,  1.76,  2.05,  2.31,  2.5, 
3.38,  3.8,  6.5,  8.0,  and  9.0  are  investigated  with  an  upstream  temperature  and  density 
of  300  K  and  6.63  xlO  6  kg/m3,  respectively,  a  Prandtl  number  of  2/3,  a  molecular 
mass  of  39.948  amu,  and  a  molecular  diameter  of  4.17  x  10-10  m  [6]. 2  The  reference 
temperature,  length,  and  time  are  Tref  =  300  K  ,  Lref  =  0.01721  m  ,  and  tref  = 
4.89  x  10-5  s  ,  respectively.  The  second-order  coupled  unsteady  BGK  and  NS  solvers 
are  used  with  the  minmod  limiter.  The  switching  parameter  value  is  Sns  =  0.001  . 
The  NS  solver  is  only  turned  on  for  the  first  7000  iterations,  at  which  point  the  NS 
and  BGK  solvers  are  coupled  together  until  the  shock  is  fully  developed  at  t  = 
0.0046  seconds  . 

The  computational  domain  is  initialized  with  the  Rankine-Hugoniot  jump  con¬ 
ditions  for  a  shockwave  [42] : 


PiUi  —  P2U2  (59) 

Pi  +  p\u\  =P2  +  P2U22  (60) 

^i  +  y  =  h2  +  Y  (61) 

where  h=  e  +  p/p  is  the  specific  enthalpy  (J/kg),  and  e  is  the  specific  internal 
energy  (J/kg).  The  computational  domain  is  split  into  two  regions:  upstream  and 
downstream  of  the  shock.  From  the  upstream  Mach  number  Mi,  density  pi,  and 
temperature  Tj,  the  rest  of  the  flow  parameters  are  determined.  The  left  and  right 
boundaries  are  set  at  their  respective  jump  conditions  for  the  duration  of  the  simula¬ 
tion.  Once  the  domain  is  initialized,  however,  the  shock  needs  to  develop,  since  it  has 
a  finite  thickness  (even  if  one  is  using  an  Euler  solver).  Therefore,  the  simulation  is 
started  with  only  the  continuum  solver  turned  on  for  a  sufficient  number  of  iterations. 
Once  the  shock  is  developed  using  the  continuum  solver,  the  Boltzmann  solver  is  also 

2The  Argon  and  Nitrogen  upstream  Mach  numbers  were  chosen  to  match  those  reported  by 
Alsmeyer  [3]. 
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turned  on  (coupled  solvers)  and  the  domain  is  allowed  to  arrive  at  a  steady-state 
solution. 

The  jump  conditions  are  derived  only  from  the  laws  of  conservation  of  mass, 
momentum,  and  energy,  respectively,  and  therefore  make  no  assumptions  about  the 
type  of  gas  being  considered  (e.g.  ideal  gas).  However,  the  Rankine-Hugoniot  jump 
conditions  do  assume  that  the  shockwave  is  a  mathematical  discontinuity  in  space 
(has  zero  thickness),  and  therefore  are  not  entirely  exact.  In  other  words,  viscosity 
and  heat  conduction  are  not  accounted  for  in  Equation  61,  which  assumes  that  the 
flow  is  adiabatic  (or  that  entropy  is  constant  along  streamlines).  In  reality,  entropy 
increases  through  the  shock  due  to  nonequilibrium,  with  the  result  that  the  entropy 
downstream  is  larger  than  the  entropy  upstream  of  the  shock.  However,  it  is  assumed 
that  these  jump  conditions  provide  adequate  boundary  conditions  for  the  flow. 

For  the  Nitrogen  cases,  Mach  numbers  Mi  =  1.53  ,  1.7,  2.0,  2.4,  2.8,  and  3.8  are 
investigated  with  an  upstream  temperature  and  density  of  300  K  and  1.0  x  10-3  kg/m3, 
respectively,  a  Prandtl  number  of  0.69,  a  molecular  mass  of  28.0  amu,  and  a  molecular 
diameter  of  4.17  x  10_1°  m  [6].  Initially,  upstream  Mach  numbers  up  to  10.0  were 
included  in  this  study,  but  due  to  the  computational  times  required  (about  2  weeks 
for  one  simulation)  as  well  as  the  memory  requirements,  it  was  determined  that  only 
Mach  numbers  up  to  3.8  would  be  considered.  The  second-order  coupled  unsteady 
3T-BGK  and  Euler  solvers  are  used  with  the  minmod  limiter.  The  NS  solver  is  not 
currently  capable  of  handling  internal  degrees  of  freedom,  so  the  Euler  solver  is  used 
instead  for  the  continuum  regions.  The  vibrational  characteristic  temperature  Ou 
is  3,371  K.3  The  rotational  collision  number  Zrot  is  5,  and  the  vibrational  collision 
number  Zm}}  is  7.9  x  107,  calculated  from  Bird’s  restatement  of  Millikan  and  White’s 

3©„  is  the  temperature  at  which  vibrational  energy  contributions  become  significant,  Zrot  is  the 
number  of  collisions  required  to  fully  transfer  rotational  energy,  and  Zvlh  is  the  number  of  molecular 
collisions  required  to  equilibrate  the  vibrational  energy. 
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where  the  constants  C-i  and  C'2  for  Ah  are  9.1  and  220.0,  respectively,  and  the  viscosity 
index  oj  is  0.74.  The  switching  parameter  value  remains  SVs  =  0.001  .  The  Euler 
solver  is  used  again  for  the  first  7000  iterations,  at  which  point  the  Euler  and  3T-BGK 
solvers  are  coupled  until  the  shock  is  fully  developed  at  t  =  0.076  seconds  . 

The  Prandtl  number  is  here  calculated  using  Eucken’s  relation  [38]: 


Pr 


4y 

7.O87  -  1.80 


(63) 


which  gives  Pr  =  0.69  for  7  =  1.40  (diatomic),  and  Pr  =  2/3  for  7  = 
1.67  (monatomic),  where  7  is  the  ratio  of  specific  heats  7  =  cp/c „  .5  These 
approximate  values  are  adequate  for  the  current  study,  even  though  it  is  recognized 
that  7  is  not  constant  through  the  temperature  ranges  considered.6 

Comparison  of  UFS  simulations  is  made  with  Alsmeyer  experimental  data  and 
with  Bird’s  educational  DSMC  code,  DSMC1S,  which  is  tailored  specifically  for  one¬ 
dimensional  shock  structure  calculations  [3] .  Alsmeyer  reported  density  data  in  Argon 
and  Nitrogen  from  a  shocktube  using  the  absorption  of  an  electron  beam,  with  an 
apparatus  similar  to  the  one  developed  by  Schmidt  [33] .  The  electron  beam  apparatus 
was  located  near  the  far  wall  of  the  shocktube,  and  is  relatively  unobtrusive.  Density, 
temperature  (translational  and  rotational),  the  axial  heat  flux  coefficient,  and  the 
partially-integrated  VDF  profiles,  shock  thickness,  and  the  density  asymmetry  factor 
are  calculated  and  presented  for  all  the  simulations  in  this  shock  structure  analysis. 


4 According  to  Jain,  Zrot  for  Aq  over  the  temperatures  considered  is  between  4  and  5  [17]. 
5Eucken  arrived  at  his  empirical  correlation  by  altering  the  result  from  kinetic  theory,  Pr  ~ 
47/(157-15)  . 

6According  to  White,  0.66  <  Pr  <  0.68  for  Argon  (a  difference  of  3%),  and  0.69  <  Pr  <  0.73  for 
N2  (a  difference  of  5%)  in  the  temperature  ranges  considered  [38]. 
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In  order  to  determine  how  well  two  density  profiles  match,  the  inverse  shock 
thickness  6,  based  on  the  maximum  density  gradient,  has  been  defined  as  [37] 


\i(— 

_  Xi  _  \dx 


xp  P2  —  pi 


(64) 


where  the  maximum  density  gradient  is  found  with  a  second-order  scheme: 


dp 

dx 


Pi+ 1  Pi— 1 
2Ax 


+  0(Ax2) 


(65) 


Bird  gives  the  upstream  equilibrium  mean  free  path  Ai  as  [6] 


_  2(5  -  2cu)(7- 2w)  /  m  \  V2  p 
1  “  15  \2^kf)  ~p 


(66) 


where  uj  is  the  viscosity  index,  p  is  the  viscosity,  p  is  the  density,  k  is  the  Boltzmann 
constant,  and  m  is  the  molecular  mass.  To  non-dimensionalize  the  experimental  data, 
Alsmeyer  calculated  Ai  based  on  the  HS  kinetic  model  (  uj  —  1/2  ): 


16  /  m  \  V2  p 
5  \2ixkT)  p 


(67) 


The  UFS  data  is  processed  using  this  formulation  of  Ai.  Bird’s  code  DSMC1S  also 
uses  the  HS  model  for  the  molecules  (  u  =  0.81  for  Argon,  to  =  0.74  for  Nitrogen), 
and  actually  outputs  Ai  (0.01304  m  for  Argon,  0.01324  m  for  Nitrogen).  By  noticing 
that 


and 


alsmeyer  _  24 
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Ai  for  the  UFS  simulations  is  calculated  to  be  0.01721  m  for  Argon  and  0.01635  m  for 
Nitrogen.  In  order  to  compare  DSMC  results  with  those  from  UFS  and  the  Alsmeyer 
data,  it  is  necessary  to  use  the  same  value  of  Ai.  Therefore,  DSMC  results  are  pre¬ 
sented  with  Ai  =  A1  aismeyer  instead  of  A1  pjrc[. 

From  a  density  profile,  one  can  construct  xp  by  extending  the  tangent  to  the 
curve  at  the  maximum  gradient  location  (  pn=0.5  )  to  the  horizontal  lines  pn  = 
1  and  pn  =  0  ,  where  the  normalized  density  pn  is 


Pn(x) 


p{x)  ~  Pi 
P2  ~  Pi 


(70) 


The  intersections  of  these  lines  are  taken  as  the  boundaries  of  the  shock,  and  so  the 
distance  between  them  on  the  x-axis  is  the  shock  thickness  {1/5),  shown  in  Figure  7. 
This  inverse  shock  thickness  does  not  take  into  account  the  skewness  of  the  density 


Figure  7:  Construction  of  shock  thickness  1/5. 
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profile,  or  any  other  behavior  upstream  or  downstream  of  the  shock  midpoint,  which 
limits  its  usefulness. 

Another  value  which  can  be  computed  is  the  density  asymmetry  factor  Q: 


=  jZ^Mx)dx 

~  fo°°  t1  “  Pn(x)\dx 


(71) 


which  is  calculated  using  a  numerical  trapezoidal  method.  Q  takes  into  account  the 
skewness  of  the  density  profile  by  dividing  the  area  under  the  curve  up  to  the  shock 
midpoint  by  the  area  between  the  curve  and  the  line  pn  =  1  beginning  at  the  shock 
midpoint.  Figure  8  gives  a  visual  representation  of  how  Q  is  calculated  (the  area  on 
the  left  divided  by  the  area  on  the  right).  A  value  of  unity  signifies  that  the  density 
profile  is  perfectly  symmetric,  while  a  value  greater  or  less  than  unity  indicates  a 
profile  skewed  to  the  left  or  right,  respectively. 


Figure  8:  Construction  of  density  asymmetry  factor  Q. 


The  overall  flow  temperature  is  not  an  output  of  UFS,  so  it  is  calculated  from 
the  non-dimensional  pressure  and  density: 

v* 

T*  =  —  (72) 

p* 

since  p*  =  p/pref  ,  P*  =  p/pref  ,  T*  =  T/Tref  ,  p  =  pRT  ,  and  pref  =  prefRTref  . 
The  normalized  temperature  profile  Tn  is  then  calculated  from 

Tn  =  (73) 

J-2  ~  J-l 

where  the  non-dimensional  superscript  *  is  dropped,  and  T\  and  T2  are  the  upstream 
and  downstream  equilibrium  temperatures,  respectively.  The  rotational  temperature 
Trot  is  an  output  of  UFS,  so  the  normalized  rotational  temperature  is  calculated  in  the 
same  way  as  the  normalized  overall  temperature  in  Equation  73.  UFS  also  outputs  the 
vibrational  temperature  Tvib,  but  DSMC1S  does  not.  Therefore,  Tvib  is  not  presented 
in  this  study. 

The  heat  flux  qx,  which  is  the  flow  of  energy  per  unit  area  per  unit  time  (W/m2), 
is  represented  using  the  heat  flux  coefficient  Cqx : 


Cqx  =  T 


~PlUlc 


(74) 


The  VDFs  of  this  study  are  two-dimensional,  but  the  VDFs  from  UFS  are  in 
terms  of  the  u-  and  u-directions,  while  the  VDFs  from  DSMC  are  in  terms  of  the  u-  and 
radial-directions.7  Therefore,  the  VDFs  from  UFS  and  DSMC  cannot  be  compared 
directly.  Instead,  the  data  is  presented  in  terms  of  a  partially-integrated  VDF  fx  [6]: 

fx  =  J  fdvdw  (75) 

7DSMC1S  assumes  that  the  non-axial  velocity  is  axially  symmetric. 
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A  few  issues  concerning  the  VDF  were  encountered  with  UFS.  First,  the  VDF 
output  from  UFS  is  not  normalized.  To  overcome  this,  the  area  under  the  curve  is 
calculated,  and  then  each  point  is  divided  by  that  total  area.  Second,  since  UFS  is  a 
relatively  new  code,  when  multiple  processors  are  used  in  parallel,  the  only  cells  that 
output  the  VDF  are  those  that  are  computed  using  the  first  processor.  This  drawback 
limits  the  results  presented  in  Chapter  IV  because  they  were  obtained  on  multiple 
processors.  Finally,  cells  flagged  as  continuum  cells  do  not  calculate  the  VDF.  So, 
when  one  is  using  UFS  to  couple  continuum  and  Boltzmann  solvers,  some  cells  of 
interest  may  not  output  a  VDF. 


3.2  Entropy- Shock  Interaction 

Two  models  for  the  entropy  spot  were  considered:  a  constant  spot,  and  a  Gaus¬ 
sian  spot.  The  constant  spot  assumes  a  discontinuity  at  the  spot  boundary  with 
constant  entropy  within  the  spot.  Even  though  a  constant  spot  may  be  easier  to 
apply  numerically,  a  more  realistic  model  can  be  implemented.  The  Gaussian  spot 
provides  a  greater  degree  of  accuracy  by  assuming  a  Gaussian  profile  of  the  form  [16] 


■r2/  2 


(76) 


shown  in  Figure  9,  where  T'  is  the  temperature  perturbation  within  the  entropy  spot 
(K),  e  is  the  perturbation  amplitude  (chosen  to  be  0.25),  and  r  is  the  radius  from  the 
center  of  the  entropy  spot  (m),  defined  as  r2  =  (x  —  xc )2  +  (y  —  yc )2  where  xc  and 
yc  are  the  x-  and  y-  coordinates  of  the  center  of  the  entropy  spot,  respectively  (m). 
The  local  density  is  then  calculated  as8 


Pi  +  Pi  — 


Pi 

7i  +  T[ 


(77) 


8Pressure  is  constant  in  an  entropy  spot,  while  density  and  temperature  fluctuate. 
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(a)  Temperature 

Figure  9:  Entropy  spot  temperature  and  density  profiles. 


also  shown  in  Figure  9.  Note  that  even  though  the  maximum  temperature  fluctuation 
is  25%,  the  minimum  density  fluctuation  is  only  —20%. 

The  grid  set-up  is  one  box  with  dimensions  20a  x  20a,  where  a  is  the  radius  of 
the  entropy  spot  (m)  (chosen  to  be  25Ai),  shown  in  Figure  10.  The  vertical  shock 
is  located  in  the  middle  of  the  domain  with  the  same  upstream  and  downstream 
boundary  conditions  as  in  the  shock  structure  study.  The  top  boundary  is  a  freestream 
boundary,  while  the  bottom  boundary  is  a  symmetry  boundary,  effectively  creating 
a  domain  that  is  20a  x  40a.  The  upstream  Mach  number  is  Mi  =  2.0  with  the 
flow  propagating  from  the  left  to  the  right.  The  center  of  the  entropy  spot  is  initially 
located  five  radii  upstream  of  the  shock  on  the  symmetry  centerline. 

The  flowheld  is  initialized  in  Argon  with  upstream  and  downstream  values  based 
on  the  Rankine-Hugoniot  shock  jump  conditions,  and  then  the  shock  is  allowed  to 
develop  with  just  the  NS  solver,  after  which  the  shock  is  allowed  to  develop  with 
the  coupled  BGK  &  NS  solvers  (  Sp  =  0.01  )  until  it  achieves  a  steady  state.9 

Once  the  flowheld  is  thus  computed  (taken  as  t  =  0  ),  the  entropy  spot  is  inserted 
upstream  of  the  shock  and  allowed  to  freely  convect  through  the  shock  using  the 

9  Pr  =  2/3  ,  molecular  mass  is  39.948  amu,  molecular  diameter  is  4.17  x  10_1°  m,  Ti  =  300 
K  ,  pi  =  6.63  x  10-6  kg/m3  ,  Ai  =  0.01721  m  . 
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Figure  10:  Grid  set-up  for  entropy-shock  simulation  at  t  =  0  (adapted  from  [16]). 

coupled  BGK  &  NS  solvers.  Figure  11  shows  the  computational  domain  (colored  by 
density)  for  Kn  =  0.1  immediately  after  the  entropy  spot  is  introduced. 

Three  Knudsen  numbers  are  investigated  in  order  to  study  the  effects  of  rarefac¬ 
tion  on  entropy-shock  interactions.  The  Knudsen  numbers  were  chosen  to  be  within 
the  continuum  (  Kn  =  0.01  ),  near-continuum  (  Kn  =  0.1  ),  and  rarefied  regimes 
(  Kn  =  1.0  ). 

Grid  studies  in  physical  space  were  conducted  on  coarse,  medium,  and  fine  grids. 
The  medium  grid  is  refined  up  to  Level  10  (  A0  =  500Ai  ,  Ai0  =  Ai/2  )  based  on 
the  spatial  gradient  of  (hip  +  In  p) .  For  the  ^-direction,  if  any  cell  has  a  gradient  of 
the  form 


d_ 

dx 


(In  p  +  lnp) 


1  dp  1  dp 
p  dx  pdx 


(78) 
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(a)  Grid  (  pmin  —  0.80  ,  Pmax  —  1.00  ) 


(b)  No  grid  (  Pmin  —  0.80  ,  Pmax  —  2.51  ) 

Figure  11:  Initialization  of  the  computational  domain,  colored  by  density  (  Kn  = 
0.1  ,  medium  physical  and  velocity  grids).  Only  a  portion  of  the  domain 
is  shown. 

larger  than  0.01  (or  of  1%),  then  it  is  refined  either  until  the  gradient  is  less  than  0.01 
or  it  is  refined  to  Level  10.  The  coarse  grid  is  refined  up  to  Level  9  (  Ag  =  Ai  ), 
and  the  fine  grid  is  refined  to  Level  11  (  An  =  Ai/4  ),  with  the  same  criterion 
on  the  gradient.  The  rarefied  simulations  (  Kn  =  1.0  )  had  to  be  computed  using 
only  the  coarse  grid.  It  is  possible  that  the  medium  and  fine  grids  were  unsuitable 
for  the  high  Knudsen  number  due  to  the  low  number  of  molecules  in  each  cell.10 
For  Kn  =  0.01  ,  the  density  profiles  agree  within  3%  between  the  fine  and  medium 
grids.  For  Kn  =  0.1  ,  the  density  profiles  agree  within  1%  between  the  coarse  and 
medium  physical  grids.  No  grid  study  was  performed  for  Kn  =  1.0  since  the  coarse 
grid  was  the  finest  grid  that  could  be  used. 

10Molecular  number  density  (number  of  molecules  per  unit  volume)  decreases  with  increasing  Kn. 
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Velocity-space  grid  studies  were  also  performed  on  coarse,  medium,  and  fine 
grids.  All  of  the  grids  have  ranges  of  — 4V%/  to  6 Vref  in  the  ^-direction,  and  —5 Vref  to 
5 Vref  in  the  u-direction.  The  coarse  grid  has  10  x  10  nodes,  the  medium  grid  has  20  x  20 
nodes,  and  the  fine  grid  has  40  x  40  nodes.  The  node  spacing  is  then  Vref,  0.5 Vref, 
and  0.251%/  for  the  coarse,  medium,  and  fine  grids,  respectively.  For  Kn  =  0.01  and 
0.1,  the  density  profiles  agree  within  0.1%  between  the  fine  and  medium  grids.  No 
velocity  grid  study  was  performed  for  Kn  =  1.0  due  to  non-physical  results  obtained 
on  the  coarse  and  fine  grids.  It  is  unclear  why  those  simulations  were  unstable. 

During  simulations,  it  was  found  that  numerical  errors  were  created  at  the 
interfaces  of  the  upper  and  lower  boundaries  and  the  shock.  The  upper  boundary  is 
far  enough  away  from  the  entropy  spot  as  to  not  be  a  concern.  However,  the  numerical 
errors  at  the  interface  of  the  shock  and  the  symmetry  boundary  (shown  in  Figure  12) 
make  it  difficult  to  analyze  the  data  right  at  the  symmetry  plane. 


Figure  12:  Numerical  error  introduced  at  the  interface  of  the  symmetry  boundary  and 
the  shock  in  density. 


Gradient  profiles  are  presented  for  density,  pressure,  and  temperature.  A  second- 
order  scheme  is  used  to  compute  the  gradients  as  before: 


db 

dx 


bi+ 1  —  &i_i 
2Ax 


+  0(  Ax2) 


(79) 


where  b  represents  density,  pressure,  or  temperature. 

No  VDFs  are  presented  since  the  only  cells  that  output  the  VDF  are  Boltzmann 
cells,  which  are  mainly  in  the  shock  structure. 
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IV.  Results 


4-1  Argon  Shock  Structure 

Nanbu  and  Watanabe  performed  a  study  on  the  one- dimensional  shock  structure 
of  a  monatomic  gas  in  1984  using  a  new  method  to  directly  simulate  the  Boltzmann 
equation  [29].  Ohwada  also  performed  a  numerical  study  of  the  Boltzmann  equation 
for  a  monatomic  gas  (HS  molecules)  and  even  reported  VDFs  [30].  However,  all  of 
their  results  agreed  well  with  those  from  Bird’s  method,  so  Bird’s  code  DSMC1S  is 
used  here  for  direct  comparison  with  the  current  results  for  both  monatomic  and 
diatomic  molecules. 

Normalized  density  profiles  in  Argon  for  various  upstream  Mach  numbers  have 
been  presented  by  Alsmeyer,  and  are  compared  here  with  the  results  from  UFS  (cou¬ 
pled  BGK  &  NS  solvers)  in  Figures  13  and  14  [3].  Also  presented  are  results 
from  UFS  using  the  NS  solver  in  the  entire  domain,  as  well  as  DSMC  results  from 
Bird’s  code  DSMC1S  [6].  BGK  agrees  closely  with  experiment  and  DSMC  up  to 
about  Mi  =  3.38  ,  after  which  BGK  deviates  significantly  by  predicting  a  thinner 
shock. 

The  inverse  shock  thickness  is  plotted  in  Figure  15;  larger  values  indicate  a 
thinner  shock,  and  vice  versa.  DSMC  overpredicts  shock  thickness,  although  it  is  the 
most  accurate  method  presented.  The  BGK  model  underpredicts  shock  thickness, 
deviating  from  experiment  at  about  Mi  =  2.5  ,  which  agrees  with  the  results 

presented  by  Schmidt,  thus  showing  that  the  BGK  model  in  UFS  performs  well  [33]. 
The  NS  solver  predicts  the  thinnest  shocks,  as  expected,  since  NS  only  allows  for 
small  deviations  from  equilibrium. 

The  density  asymmetry  factor  is  plotted  in  Figure  16  with  experimental  values 
from  Alsmeyer  and  Schmidt  [3,33].  The  experimental  results  indicate  that  density 
profiles  skew  more  to  the  left  with  increasing  Mi.  DSMC  follows  experiment  fairly 
well,  while  the  coupled  UFS  solver  (BGK  &  NS)  profiles  are  skewed  much  more  to 
the  left.  The  NS  solver  shows  a  relatively  constant  value  for  Q  with  the  profiles  being 
skewed  to  the  left  for  all  upstream  Mach  numbers. 
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(c)  Ml  =  2.05 


(d)  Mi  =  2.31 


Figure  13:  Density  profiles  in  Argon.  BGK  &  NS  (blue);  NS  (green);  DSMC  (red); 
Alsmeyer  (black). 


The  translational  temperature  Ttr  profiles  for  Argon  are  shown  in  Figures  17  and 
18  for  the  UFS  coupled  solver  (BGK  &  NS),  the  UFS  NS  solver,  and  DSMC.  Since 
Argon  is  monatomic,  it  has  no  internal  energy  modes  (besides  electronic,  which  is  ne¬ 
glected  in  this  study) ,  and  therefore  the  only  molecular  energy  that  can  be  transferred 
is  kinetic  (or  translational),  which  means  that  these  profiles  also  represent  the  overall 
temperatures  in  the  flow.  The  temperature  profiles  have  midpoint  values  upstream 
of  the  midpoint  of  the  shock  based  on  density.  The  NS  profiles  show  larger  gradients 
with  increasing  M\  and  provide  limiting  cases  for  comparison.  Up  to  Mi  =  3.8  , 


46 


(a)  Mi  =  3.38 


-8  -6  -4  -2  0  2  4  6  8 


xA-1 

(b)  Mi  =  3.8 


(c)  Mi  =  6.5 


(d)  Mi  =  9.0 


Figure  14:  Density  profiles  in  Argon.  BGK  &  NS  (blue);  NS  (green);  DSMC  (red); 
Alsmeyer  (black). 


DSMC  and  BGK  &  NS  show  good  agreement,  after  which  DSMC  predicts  a  higher 
relaxation  time. 

The  heat  flux  coefficient  C(Jx  is  shown  in  Figures  19  and  20.  The  minima  of  the 
coupled  UFS  solver  (BGK  &  NS)  and  the  DSMC  code  are  within  2%  for  all  values 
of  Mi,  with  the  profiles  matching  well  up  to  M\  =  3.8  .  At  higher  upstream  Mach 
numbers,  DSMC  shows  higher  relaxation  times,  an  earlier  onset  of  nonequilibrium, 
and  minima  which  occur  upstream  of  the  BGK  &  NS  minima. 

The  partially-integrated  VDFs  (/.,,)  are  given  in  Figures  21  and  22  for  the  DSMC 
and  BGK  &  NS  solvers.  Up  to  Mi  =  2.5  ,  the  profiles  show  very  good  agreement. 
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Figure  15:  Inverse  shock  thickness  S  for  Argon. 


At  Mi  =  3.38  and  3.8,  the  profiles  begin  to  show  some  deviation  from  one  another 
with  the  BGK  &  NS  profiles  being  skewed  slightly  more  to  the  right  than  the  DSMC 
profiles. 

Even  greater  deviation  is  exhibited  for  Mi  =  6.5  ,8.0,  and  9.0.  This  can  be 
explained  by  a  phenomenon  encountered  during  simulation,  namely  that  the  shock  in 
UFS  began  to  move  upstream.1  The  other  results  (i.e.  density,  temperature,  and  heat 
flux)  are  unaffected  by  a  moving  shock  since  they  all  reference  the  density  midpoint 
of  the  shock.  However,  the  VDFs  are  output  from  UFS  at  specific  locations  in  the 
domain,  so  the  errors  seen  at  Mi  =  6.5  ,  8.0,  and  9.0  are  apparent  since  the  midpoint 
of  the  shock  is  not  located  at  the  center  of  the  cell  (but  rather  between  two  cells). 

1 A  moving  shock  was  not  observed  for  the  other  Mach  numbers. 


48 


Figure  16:  Density  asymmetry  Q  in  Argon.  The  solid  curve  is  from  Alsmeyer  [3]. 

From  the  previous  results,  it  is  clear  that  for  these  Mach  numbers,  DSMC  predicts 
greater  nonequilibrium  farther  upstream  than  BGK  &  NS,  indicating  that  the  DSMC 
VDF  evolves  more  slowly  than  the  BGK  &  NS  VDF  as  the  gas  convects  through  the 
shock.  Therefore,  the  BGK  &  NS  VDF  at  x/Xi  =  — 1.5  has  a  higher  maximum  than 
the  DSMC  VDF,  and  a  lower  maximum  than  the  DSMC  VDF  at  x/Xi  =  0.0  . 

The  total  number  of  hours  required  to  run  each  simulation  (including  the  grid 
independence  studies)  is  shown  in  Figure  23  versus  the  upstream  Mach  number  Mi. 
At  Mi  =  1.55  the  time  step  is  0.229  /rs,  while  at  Mi  =  9.0  the  time  step  is  0.161  /rs. 

For  UFS,  the  run  times  are  relatively  constant  up  to  Mi  =  6.5  ,  at  which  point 
they  increase  exponentially,  as  shown  on  the  logarithmic  scale.  The  DSMC  simulations 
took  from  1  to  2  hours  to  complete,  while  the  UFS  simulations  required  up  to  122 
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=  (T-T 1 )/ (T  2-T 1 )  T  =  (T-T^T^)  T  =  (T-T^T^) 


1.2 


1.2 


(a)  Mi  =  1.55 


(b)  M\  =  1.75 


xA1  x/k 

(c)  Mi  =  1.76  (d)  Mi  =  2.05 


(e)  Mi  =  2.31  (f)  Mi  =  2.5 

Figure  17:  Translational  temperature  profiles  in  Argon.  BGK  &  NS  (blue 
NS  (green);  DSMC  (red). 
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=  (T-T1)/(T2-T.|)  Tn  =  (T_Ti  )/(T2_Ti ) 


(a)  Mi  =  3.38  (b)  M1  =  3.8 


(c)  Mi  =  6.5  (d)  Mi  =  8.0 


(e)  Mi  =  9.0 


Figure  18:  Translational  temperature  profiles  in  Argon.  BGK  &  NS  (blue); 
NS  (green);  DSMC  (red). 
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(a)  Mi  =  1.55 


(b)  Mi  =  1.75 


(c)  Mi  =  1.76 


(d)  Mi  =  2.05 


(e)  Mi  =  2.31 


(f)  Mi  =  2.5 


Figure  19:  Heat  flux  coefficient  profiles  for  qx  in  Argon.  BGK  &  NS  (blue);  NS  (green) 
DSMC  (red). 
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(a)  Mi  =  3.38  (b)  M1  =  3.8 


(c)  Mi  =  6.5 


(d)  Mi  =  8.0 


(e)  Mi  =  9.0 

Figure  20:  Heat  flux  coefficient  profiles  for  qx  in  Argon.  BGK  &  NS  (blue);  NS  (green) 
DSMC  (red). 
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(a)  Mi  =  1.55 


(c)  Mi  =  1.76 


(e)  Mi  =  2.31 

Figure  21:  Partially- integrated  VDFs  ir 
—  13.5  (blue);  x/\\  =  —1.5 


(b)  Mi  =  1.75 


(d)  Mi  =  2.05 


(f)  Mi  =  2.5 


Argon.  -  BGK  &  NS;  — ,  DSMC.  x/X1 
(green);  x/\\  =  0.0  (red) 
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(c)  Ml  =  6.5  (d)  Ml  =  8.0 


(e)  Mi  =  9.0 

Figure  22:  Partially-integrated  VDFs  in  Argon.  — ,  BGK  &  NS;  — ,  DSMC.  x/\\ 
—  13.5  (blue);  x/\\  =  —1.5  (green);  x/X\  =  0.0  (red) 


55 


1000 


1 


9 


-•-Coarse  phys  grid 

♦  Medium  phys  grid 

♦  Fine  phys  grid 
— Coarse  vel  grid 

♦  Fine  vel  grid 

♦  BGKonly 

♦  DSMC 


Figure  23:  Total  run  times  (in  hours)  versus  Mi  for  the  Argon  simulations  on  a  semi¬ 
log  scale. 


hours  (61  x  longer)  for  the  medium  grid.  It  must  be  noted  that  the  DSMC  solver 
uses  a  one- dimensional  grid,  whereas  the  UFS  simulations  require  a  two-dimensional 
grid.  Also,  the  code  DSMC1S  was  specifically  designed  to  compute  one-dimensional 
shock  structures,  while  UFS  is  a  general  code  meant  to  compute  flowhelds  much  more 
complex  than  just  a  one- dimensional  shock.  Therefore,  it  makes  sense  that  UFS 
requires  more  computationally  expensive  than  DSMC1S. 


4-2  Nitrogen  Shock  Structure 

Normalized  density  profiles  in  Nitrogen  for  various  upstream  Mach  numbers 
have  been  presented  by  Alsmeyer,  and  are  compared  here  with  the  results  from  UFS 
(coupled  3T-BGK  &  Euler  solvers)  in  Figure  24  [3].  Also  presented  are  results  from 
UFS  using  the  Euler  solver  in  the  entire  domain,  as  well  as  DSMC  results  from  Bird’s 
code  DSMC1S  [6].  One  immediately  notices  that  the  Euler  results  are  atypical.  A 
regular  Euler  solver  predicts  a  very  thin  shock,  whereas  UFS’s  Euler  solver  shows  a 
relaxation  region  downstream  of  the  midpoint  of  the  shock.  This  discrepancy  may 
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(ld-2d)/(Ld  -  d)  =  ud  (ld-2d)/(Ld  -  d)  =  ud  (ld-2d)/(ld  -  d)  = 


(a)  Mi  =  1.53 


(b)  Mi  =  1.7 


(c)  Mi  =  2.0 


(d)  Mi  =  2.4 


(e)  Mi  =  2.8 


(f)  Mi  =  3.8 


Figure  24:  Density  profiles  in  Nitrogen.  3T-BGK  &  Euler  (blue);  Euler  (green); 
DSMC  (red);  Alsmeyer  (black). 
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be  explained  by  the  fact  that  the  Euler  solver  in  UFS  is  capable  of  accounting  for 
internal  energies,  while  a  classical  Euler  solver  is  not.  The  developers  have  not  yet 
published  how  the  Euler  solver  was  modified,  so  further  explanation  of  these  results 
must  be  deferred.  The  DSMC  data  closely  match  the  experimental  results,  and  the 
coupled  solver  predicts  thicker  shocks  and,  consequently,  higher  relaxation  times. 
Agreement  with  the  coupled  solver  and  experiment  improve  with  increasing  upstream 
Mach  number. 

The  inverse  shock  thickness  is  plotted  in  Figure  25;  larger  values  indicate  a 
thinner  shock,  and  vice  versa.  DSMC  agrees  well  with  Alsmeyer’s  curve,  even  though 


Figure  25:  Inverse  shock  thickness  8  for  Nitrogen. 


it  overpredicts  shock  thickness  for  Mi  >  4.0  .  The  3T-BGK  &  Euler  simulations 
show  a  similar  trend  as  the  BGK  &  NS  simulations  in  Figure  15,  since  they  predict 
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thinner  shocks  than  experiment  for  Mi  >  3  .  The  Euler  results  are  not  shown  due 
to  the  atypical  shock  profiles  given  in  Figure  24. 

The  density  asymmetry  factor  is  plotted  in  Figure  26  with  experimental  values 
from  Alsmeyer  for  reference  [3].  Even  though  Nitrogen  simulations  are  here  being 


Figure  26:  Density  asymmetry  Q  in  Nitrogen. 


compared  with  Argon  simulations,  DSMC  closely  follows  the  Alsmeyer  curve,  while 
the  coupled  UFS  solver  (3T-BGK  &  Euler)  agrees  with  DSMC  for  Mi  >  2.0  . 

As  Nitrogen  flows  through  the  shock,  the  gas  experiences  translational,  rota¬ 
tional,  and  vibrational  nonequilibrium.  The  translational  energy  of  a  molecule  is 
completely  equilibrated  after  only  one  or  two  collisions,  while  it  takes  five  and  79 
million  collisions  for  the  rotational  and  vibrational  energies  to  be  completely  trans¬ 
ferred,  respectively.  As  a  result,  Ttr  not  only  equilibrates  before  Trot,  but  Ttr  may 
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be  temporarily  higher  than  its  equilibrium  value.  This  phenomenon  occurs  because 
it  takes  many  more  collisions  for  energy  to  be  transferred  to  the  rotational  and  vi¬ 
brational  modes,  and  yet  the  flow  energy  must  go  somewhere.  Hence,  the  energy 
that  will  eventually  go  into  rotational  and  vibrational  excitation  is  temporarily  stored 
as  translational  energy.  The  translational  and  rotational  temperature  profiles  are 
shown  in  Figure  27.  The  coupled  3T-BGK  &  Euler  solvers  predict  higher  relaxation 
times  for  the  rotational  energy  than  does  the  DSMC  solver.  As  a  result,  the  coupled 
translational  temperatures  are  higher  than  those  for  DSMC  by  as  much  as  2%. 2 

Heat  flux  coefficient  C<h.  profiles  are  shown  in  Figure  28.  The  Euler  simulations 
are  not  given  since  qx  is  everywhere  zero.  For  the  upstream  Mach  numbers  presented, 
the  DSMC  and  3T-BGK  &  Euler  solvers  differ  at  their  minima  by  as  much  as  15%. 
The  3T-BGK  &  Euler  solver  also  experiences  a  maximum  between  1  <  x/\\  <  5  , 
which  is  not  apparent  from  the  DSMC  data.  It  is  unclear  whether  this  result  is  a 
numerical  artifact  of  the  3T-BGK  kinetic  scheme,  or  whether  DSMC  is  incorrectly 
predicting  the  heat  flux,  especially  since  there  are  arguments  for  both.3 

The  partially-integrated  VDFs  (fx)  are  given  in  Figure  29.  The  coupled  3T- 
BGK  &  Euler  solver  shows  very  good  agreement  with  the  DSMC  code.  However,  it 
is  expected  that  a  more  refined  velocity  grid  would  give  even  better  agreement,  since 
the  grid  spacing  here  is  0.5V%/. 

The  total  number  of  hours  required  to  run  each  simulation  (including  the  grid 
independence  studies)  is  shown  in  Figure  30  versus  the  upstream  Mach  number  M\. 
At  Mi  =  1.53  the  time  step  is  5.05  /ts,  while  at  Mi  =  3.8  the  time  step  is  2.81  /ts. 

For  UFS,  the  run  times  are  relatively  constant  up  to  Mi  =  2.0  ,  at  which 
point  they  increase  exponentially,  as  shown  on  the  logarithmic  scale.  The  DSMC 

2  Comparison  is  made  between  dimensional  temperatures. 

3In  Figure  27,  the  translational  temperature  experiences  a  maximum  above  the  downstream 
equilibrium  temperature,  thus  requiring  a  positive  heat  flux,  justifying  the  3T-BGK  &  Euler  results. 
Alternatively,  the  energy  captured  in  the  translational  temperature  when  it  peaks  goes  into  the 
rotational  energy  modes  as  the  rotational  temperature  equilibrates,  which  requires  only  the  heat 
flux  predicted  by  the  DSMC  results. 
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(a)  Mi  =  1.53 


(b)  Mi  =  1.7 
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(c)  Mi  =  2.0 


(d)  Mi  =  2.4 
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(e)  Mi  =  2.8 


(f)  Mi  =  3.8 


Figure  27:  Translational  and  rotational  temperature  profiles  in  Nitrogen.  — ,  Ttr;  - 
Trot ;  3T-BGK  &  Euler  (blue);  DSMC  (red). 
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(a)  Mi  =  1.53  (b)  Mi  =  1.7 


(c)  Mi  =  2.0 


(d)  Mi  =  2.4 


(e)  Mi  =  2.8  (f)  Mi  =  3.8 


Figure  28:  Heat  flux  coefficient  profiles  for  qx  in  Nitrogen.  3T-BGK  &  Euler  (blue) 
DSMC  (red). 
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(a)  Mi  =  1.53  (b)  Mi  =  1.7 


(c)  Mi  =  2.0  (d)  Mi  =  2.4 


(e)  Mi  =  2.8  (f)  Mi  =  3.8 


Figure  29:  Partially-integrated  VDFs  in  Nitrogen. 

DSMC.  x/\\  =  —13.5  (blue);  x/\\  =  —1.5 


3T-BGK  &  Euler;  — , 
(green);  x/\\—  0.0  (red) 
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Figure  30:  Total  run  times  (in  hours)  versus  M\  for  the  Nitrogen  simulations  on  a 
semi-log  scale. 


simulations  took  from  1  to  2  hours  to  complete,  while  the  UFS  simulations  required 
up  to  42  hours  (24 x  longer)  for  the  medium  grid.  The  coarse  and  fine  velocity  grids 
bracket  the  UFS  simulation  times,  showing  that  the  most  important  factor  in  run 
times  is  the  number  of  nodes  in  velocity  space. 


4-  3  Entropy- Shock  Interaction 

The  midpoint  of  the  shock  based  on  density  moves  as  the  entropy  spot  con- 
vects  through  the  shock,  and  the  entire  shock  actually  is  bowed,  however  slightly. 
Figure  31  gives  the  change  in  shock  location  for  all  three  Knudsen  numbers,  where 
the  non-dimensional  time  r  equals  zero  when  the  center  of  the  entropy  spot  passes 
through  x/\i  —  0  .  Since  the  center  of  the  entropy  spot  is  5a  upstream  of  the  shock, 
r  is  calculated  as  [12] 

t-  5  — 

7  =  ’  <80) 
Ui 
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(a)  All  Kn  (b)  Kn  =  0.01 


(c)  Kn  =  0.1  (d)  Kn  =  1.0 

Figure  31:  Shock  locations  as  the  simulation  progresses  in  time  (r)  for  all  three  Knud- 
sen  numbers. 

where  t  is  the  time  (s),  a  is  the  radius  of  the  entropy  spot,  and  u\  is  the  upstream 
^-direction  velocity  (m/s).  Qualitatively,  the  shock  location  profiles  exhibit  the  same 
behaviors,  just  on  different  scales.  Even  at  Kn  =  0.01  ,  which  is  well  into  the 

continuum  regime,  the  shock  experiences  a  restoring  force  towards  x/\\  =  0  ,  albeit 
a  small  one.  It  is  unclear  why  the  shock  does  not  return  to  its  equilibrium  location  (as 
a  completely  normal  shock)  at  Kn  =  0.01  .A  possible  explanation  may  be  that  the 
pressure  wave,  which  emanates  from  the  entropy  spot  after  it  has  impinged  with  the 
shock,  is  more  powerful  at  lower  Knudsen  numbers,  causing  the  shock  to  be  bowed 
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upstream.  If  this  is  the  case,  then  the  pressure  wave  at  Kn  =  0.1  and  1.0  is  much 
weaker  so  as  to  allow  the  shock  to  almost  completely  be  restored  to  a  normal  shock. 

The  density,  pressure,  and  temperature  profiles,  as  well  as  their  gradient  profiles, 
are  shown  in  Figures  33,  34,  and  35  for  three  different  Knudsen  numbers,  Kn  =  0.01  , 
0.1,  and  1.0,  respectively. 

For  Kn  =  0.01  ,  the  density  fluctuation  is  initially  unaffected  by  the  shock,  with 
a  post-shock  strength  of  —20%  at  r  =  1.0  ,  damping  out  to  —14%  by  r  =  20.3  .  It  is 
possible  that  with  such  a  low  Knudsen  number,  the  shock  thickness  is  not  allowing  the 
density  fluctuation  sufficient  time  to  change.  The  temperature  fluctuation,  however  is 
initially  affected  by  the  shock,  which  has  a  post-shock  strength  of  only  20%  (compared 
to  a  pre-shock  strength  of  25%)  at  r  =  1.0  ,  damping  out  to  14%  at  r  =  20.3  .  It 
is  interesting  to  note  that  even  though  the  entropy  spot  does  not  induce  any  pressure 
fluctuations  upstream  of  the  shock,  by  convecting  through  the  shock,  the  interaction 
creates  a  circular  pressure  wave  (see  [12])  centered  about  the  symmetry  plane.  The 
pressure  profile  in  Figure  33  shows  a  minimum  of  —5%  at  r  =  1.0  which  dampens 
to  a  minimum  of  —2%  at  r  =  5.8  .  For  r  >  5.8  ,  the  pressure  profile  shows  that 
there  is  a  rarefaction  region  at  the  symmetry  line.  The  numerical  error  introduced 
by  the  shock-symmetry  boundary  interface  is  apparent  here  near  x/a  =  0  since  the 
flow  parameters  are  relaxing  from  values  above  the  downstream  equilibrium  values. 

For  Kn  =  0.1  ,  the  shock  thickness  is  larger  than  the  simulation  with  Kn  = 
0.01  .  As  a  result,  the  density  fluctuation  is  more  affected  initially  by  the  shock,  with 
a  post-shock  strength  of  —12%  at  r  =  0.8  ,  damping  out  to  —7%  by  r  =  18.3  . 
The  strength  of  the  temperature  fluctuation  is  cut  in  half  to  12%  at  r  =  5.2  , 

damping  out  to  8%  at  r  =  20.3  .  The  pressure  shows  a  minimum  of  —5%  at  r  = 
0.8  which  dampens  to  a  minimum  of  —1.5%  at  r  =  5.2  .  For  r  >  5.2  ,  the 
pressure  profile  indicates  that  there  is  actually  a  compression  region  after  the  shock 
at  Kn  =  0.1  ,  which  is  in  contrast  to  the  rarefaction  region  predicted  for  Kn  =  0.01  . 
The  numerical  error  introduced  by  the  shock-symmetry  boundary  interface  is  again 
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apparent  here  near  x/a  =  0  since  the  flow  parameters  are  relaxing  from  values  above 
the  downstream  equilibrium  values. 

For  Kn  =  1.0  ,  the  shock  thickness  is  larger  still,  apparent  from  the  attenuation 
of  the  density  fluctuation  downstream  of  the  shock  (—7%  at  r  —  2.1  ,  dampened 
to  —1.4%  at  r  =  9.6  ).  The  numerical  error  from  the  shock-symmetry  boundary  is 
lessened  at  this  Knudsen  number  (possibly  due  to  the  thicker  shock),  and  so  the  flow 
parameters  (except  for  pressure)  are  qualitatively  showing  correct  behavior  directly 
downstream  of  the  shock  (relaxing  from  values  below  the  downstream  equilibrium 
values).  The  strength  of  the  temperature  fluctuation  is  now  only  4%  at  r  =  0.8  , 
damping  out  to  2%  by  r  =  9.6  .  The  pressure  actually  is  exhibiting  values  above  the 
downstream  equilibrium  value  by  3%  at  r  =  0.8  ,  and  4%  by  r  =  9.6  ,  suggesting 
a  compression  region  downstream  of  the  shock. 

The  total  number  of  hours  required  to  run  each  simulation  (including  the  grid 
independence  studies)  is  shown  in  Figure  32  versus  the  Knudsen  number  Kn. 


-♦-Coarse  phys  grid 
-■-Medium  phys  grid 
-a- Fine  phys  grid 
^Coarse  vel  grid 
^Fine  vel  grid 


Kn 


Figure  32:  Total  run  times  (in  hours)  versus  Kn.  The  coarse  and  fine  velocity  grids 
use  a  medium  physical  grid,  and  the  coarse  and  fine  physical  grids  use  a 
medium  velocity  grid. 
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(a)  Density 


(b)  Density  gradients 


(c)  Temperature  (d)  Temperature  gradients 


(e)  Pressure 


(f)  Pressure  gradients 


Figure  33:  Profiles  for  Kn  =  0.01 
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(a)  Density 


(b)  Density  gradients 


(c)  Temperature  (d)  Temperature  gradients 


(e)  Pressure 


(f)  Pressure  gradients 


Figure  34:  Profiles  for  Kn  =  0.1 
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(c)  Temperature 


(d)  Temperature  gradients 


(e)  Pressure  (f)  Pressure  gradients 


Figure  35:  Profiles  for  Kn  =  1.0 
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V.  Conclusions 


The  schemes  implemented  in  UFS  (Euler,  NS,  BGK,  and  3T-BGK)  provide  accu¬ 
rate  results  for  one-dimensional  simulations  of  shock  structure  in  both  Argon  and 
diatomic  Nitrogen.  Profiles  of  density,  temperatures,  heat  flux  coefficient,  and  the 
partially-integrated  VDF,  as  well  as  the  shock  thickness  and  the  density  asymmetry 
factor,  show  good  agreement  when  compared  with  experimental  results  reported  by 
Alsmeyer  and  with  numerical  simulations  on  Bird’s  DSMC1S  code  [3].  The  Euler 
scheme  in  UFS  gives  atypical  results,  possibly  due  to  its  inclusion  of  internal  energies. 
Coupling  of  continuum  and  kinetic  solvers  allows  one  to  simulate  a  large  flowfield 
while  keeping  computational  cost  down  in  flows  with  both  continuum  and  rarefied 
regions.  Further  study  of  shock  structure  may  include  the  other  models  and  solution 
methods  included  in  UFS  (e.g.  DNS  of  the  collision  integral,  and  various  intermolec- 
ular  potential  models),  and  should  utilize  a  one-dimensional  grid  in  order  to  further 
reduce  the  computational  cost  required.  Studies  of  additional  continuum  breakdown 
parameters  may  be  made  with  comparison  to  this  work. 

Entropy-shock  interactions  continue  to  be  an  area  of  research,  with  the  work 
presented  here  extending  the  analysis  to  NS  and  kinetic  schemes.  The  shock  moves 
upstream  due  to  entropy-shock  interactions,  and  seems  to  recover  to  a  normal  shock 
for  higher  Knudsen  numbers  (  Kn  =  0.1  and  1.0).  Also,  for  Kn  =  0.01  ,  a 

rarefaction  region  is  observed  downstream  of  the  shock,  while  for  Kn  =  0.1  and 
1.0  a  compression  region  is  given.  Numerical  errors  introduced  where  the  shock  and 
the  symmetry  boundary  interface  create  uncertainty  in  the  results  presented.  Further 
study  should  investigate  the  cause  of  this  effect,  and  either  use  a  full  domain  (instead 
of  half  of  the  domain  with  a  symmetry  boundary)  or  modify  the  symmetry  boundary 
condition  to  eliminate  the  errors.  As  CFDRC  improves  the  capabilities  of  UFS,  es¬ 
pecially  the  visualization  software,  the  pressure  waves  downstream  of  the  shock  may 
be  studied  using  pressure  contours.  VDFs  may  be  studied  by  using  a  smaller  value 
for  the  continuum  breakdown  parameter,  or  by  performing  the  simulations  with  an 
uncoupled  kinetic  solver. 
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Appendix  A:  List  of  Symbols 


Symbol  Page 

a  radius  of  the  entropy  spot  (m)  41 

A\r  local  constants  (  i  —  1,  2,  3,  4,  5  )  25 

A1  local  constants  (  *  =  1,  2,  3,  4,  5  )  25 

a\  local  constants  (  i  —  1,  2,  3,  4,  5  )  25 

a\  local  constants  (  i  —  1,  2,  3,  4,  5  )  25 

Cj  molecular  velocity  (m/s) .  12 

Ci  average  molecular  velocity,  or  macroscopic  gas  velocity  (m/s)  .  12 

c'  post-collisional  molecular  velocity  (m/s) .  17 

cp  specific  heat  at  constant  pressure  (J/kg-K)  .  26 

Ci  thermal  velocity  (m/s) .  12 

CQx  heat  flux  coefficient  .  39 

d  molecular  diameter  (m) .  17 

dVc  volume  element  in  velocity  space  (m3/s3)  15 

dVx  volume  element  in  physical  space  (m3) .  15 

S  inverse  shock  thickness  .  36 

Ao  dimension  for  base  cells  (m)  (refinement  Level  0) .  29 

A2  cell  dimension  for  cells  of  refinement  Level  2  (m) .  29 

e  specific  internal  energy  ( J /kg)  .  33 

E  energy  (J)  22 

e  perturbation  amplitude .  40 

/  normalized  VDF .  11 

feg  equilibrium  (or  Maxwellian)  VDF  .  12 

/o  VDF  at  an  initial  moment  (  t  —  0  ) .  24 

fk7l  VDF  value  on  the  cell  face  .  20 

fi9  equilibrium  VDF  on  the  cell’s  left  face .  23 

//9  equilibrium  VDF  on  the  cell’s  right  face .  23 
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Symbol  Page 

fx  partially-integrated  VDF  in  the  axial  direction .  39 

Ft  molecular  acceleration  (m/s2) .  16 

F^i/2^  flux  on  the  cell  face  along  the  x-direction .  22 

g  relative  molecular  speed  (m/s) .  17 

G  /J+i/2  /.flux  on  the  cell  face  along  the  ^-direction .  22 

7  ratio  of  specific  heats  .  35 

h  non-dimensional  local  cell  size  .  25 

h  specific  enthalpy  (J/kg) .  33 

H'-'hk+1,2f{ux  on  the  cell  face  along  the  x-direction .  22 

H[£]  step  function .  23 

i  node  in  velocity  space  (subscript) .  20 

i  cell  node  in  physical  space  in  the  x-direction  (subscript)  ....  22 

I (/,  /)  collision  integral .  20 

j  cell  number  in  physical  space  (subscript) .  20 

j  cell  node  in  physical  space  in  the  ^-direction  (subscript)  ....  22 

k  Boltzmann  constant  (1.38  x  10-23  kg  •  m2/s2  •  K) .  12 

k  time  index  (superscript) .  20 

k  cell  node  in  physical  space  in  the  ^-direction  (subscript)  ....  22 

k  coefficient  of  thermal  conductivity  (W/m-K)  .  26 

Kr  rotational  degrees  of  freedom .  27 

Kv  vibrational  degrees  of  freedom  .  27 

Kn  Knudsen  number .  5 

L  characteristic  flow  length  (m) .  5 

Lref  reference  length  (m) .  29 

A  mean  free  path  (m)  4 

m  molecular  mass  (kg) .  12 

Mn  nth  moment  of  the  VDF .  11 

/j,  dynamic  viscosity  (N-s/m2) .  24 
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Symbol  Page 

n  number  density  (number  of  molecules  per  unit  volume,  1/m3)  .  14 

n  time  step  (superscript)  .  22 

n  cell  face  unit  outward  normal  vector .  20 

N  number  of  molecules .  17 

u  collision  frequency  (Hz,  or  1/s) .  18 

u  viscosity  index .  35 

p  pressure  (N/m2) .  13 

Pr  Prandtl  number .  26 

<f>  inverse  collision  integral .  21 

<f>(Ci)  velocity  component  distribution  function .  12 

•0  collision  invariants .  22 

qL  heat  flux  (W/m2)  13 

Q  any  quantity  that  is  a  function  of  velocity .  11 

Q  heat  flux  (W/m2)  26 

Q  density  asymmetry  factor .  38 

r  radius  from  the  center  of  the  entropy  spot  (m)  40 

r  position  vector  in  physical  space  (m) .  20 

R  specific  gas  constant  ( J /kg  •  K) .  14 

M3  set  of  all  real  numbers  in  three-dimensional  space .  22 

p  density  (kg/m3)  13 

pn  normalized  density .  37 

S face  face  surface  area  (m2) .  20 

t  time  (s) .  65 

At  non-dimensional  time  step .  25 

ti  dummy  integration  variable  for  time  (s) .  24 

T  temperature  (K) .  12 

T  non-dimensional  temperature .  25 

Tn  normalized  temperature .  39 
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Symbol  Page 

Tref  reference  temperature  (K) .  32 

Trot  rotational  temperature  (K) .  27 

Ttr  translational  temperature  (K) .  27 

Tv n  vibrational  temperature  (K) .  27 

T'  temperature  perturbation  (K) .  40 

Teq  equilibrium  temperature  (K)  27 

r  intercollision  relaxation  time  (s)  24 

r  non-dimensional  time .  64 

Tij  shear  stress  (N/m2)  .  13 

tu  local  BGK  relaxation  time  (s) .  18 

0,;  vibrational  characteristic  temperature  (K) .  34 

u  flow  velocity  in  the  x-direction  (m/s) .  22 

Ui  upstream  velocity  (m/s) .  65 

U  non-dimensional  flow  velocity .  25 

v  flow  velocity  in  the  ^-direction  (m/s) .  22 

V  cell  volume  (m3) .  20 

Vref  reference  velocity  (m/s) .  32 

w  flow  velocity  in  the  ^-direction  (m/s) .  22 

xc  x-coordinate  of  the  center  of  the  entropy  spot  (m)  40 

X[  trajectory  of  a  particle  (m) .  24 

£, max  non-dimensional  molecular  velocity .  26 

normal  velocity  to  the  cell  face  (m/s)  .  25 

£  velocity  vector  (m/s)  20 

velocity  space  nodes .  20 

velocity  space  cell  size  (m/s)  20 

yc  ^-coordinate  of  the  center  of  the  entropy  spot  (m)  40 

Y"ijk  cell- averaged  value  at  time  step  n  .  22 

Zrot  rotational  collision  number .  34 
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Symbol  Page 

Zmij  vibrational  collision  number .  34 

Q  pre-collisional  molecular  velocity  (m/s)  17 

Q  post-collisional  molecular  velocity  (m/s) .  17 

*  intermediate  time  level  (superscript) .  20 


76 


Appendix  B:  List  of  Abbreviations 

Abbreviation  Page 

3T-BGK  Three- Temperature  BGK  .  26 

AFRL  Air  Force  Research  Laboratory .  1 

BGK  Bhatnagar-Gross-Krook  collision  model .  18 

CFD  Computational  Fluid  Dynamics .  2 

CFDRC  CFD  Research  Corporation .  7 

CFL  Courant-Friedrichs-Lewy .  25 

DNS  Direct  Numerical  Simulation .  19 

DoD  Department  of  Defense .  1 

DSMC  Direct  Simulation  Monte  Carlo .  6 

HiFire  Hypersonic  International  Flight  Research  Experimentation  1 

HS  Hard  Sphere  .  19 

ITAR  International  Traffic  in  Arms  Regulations  .  27 

MPI  Message-Passing  Interface .  19 

NIWA  National  Institute  of  Water  and  Atmospheric  research  .  .  19 

NS  Navier-Stokes .  2 

UFS  Unified  Flow  Solver .  7 

USAF  United  States  Air  Force .  1 

VDF  Velocity  Distribution  Function  .  3 
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